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SUMMARY 


STABCAR is a computer program that can be used to determine the characteristic 
roots of flexible, actively controlled aircraft, including the effects of unsteady 
aerodynamics, A modal formulation is employed to characterize the aircraft, The 
control system representation is input as a matrix of transfer functions. Roots are 
determined through the use of determinant iteration utilizing either oscillatory or 
approximate fully unsteady aerodynamic forces. Provisions are included to allow 
automatic variation of velocity, density, altitude, or feedback gains. 

The mathematical model employed for the aircraft is described, a flowchart is 
provided, detailed description is given of the function of key program elements and 
parameters that direct program operation, program inputs are defined, and sample 
inputs and outputs are presented. 

STABCAR can be executed in either a batch or an interactive mode. It is written 
specifically for use on CDC® CYBER 175 eguipment. Modification would be required for 
operation on other equipment. 


INTRODUCTION 

The equations of motion of flexible aircraft contain aerodynamic force terms 
which, when expressed in the Laplace domain, are transcendental functions of the 
Laplace variable. Determination of the characteristic roots of the system, there- 
fore, requires the use of iterative or graphical procedures. Alternatively, one may 
approximate the aerodynamic force terms with rational functions of the Laplace vari- 
able and convert the equations of motion to a standard eigenvalue problem (refs. 1 
and 2). 

This paper describes a computer program which utilizes the method of determinant 
iteration to solve for the characteristic roots of flexible, actively controlled 
aircraft. References 3 through 6 describe other programs for iterative determination 
of characteristic roots. Reference 3, which was closely followed in the development 
of STABCAR, describes results obtained with a code developed by a major airframe 
manufacturer; reference 4 describes methods employed in Ehgland. Neither code is 
readily available and no reference is made to a user's guide for the programs. Ref- 
erences 5 and 6, however, describe NASTRAN®, a publically available program. The 
NASTRAN program was developed for batch-type operations. In contrast, STABCAR, 
described herein, is tailored for interactive execution and contains a number of 
features which facilitate the interactive study of the effects upon stability of 
variations in control-system parameters. 

STABCAR is a module of the ISAC (Interaction of Structures, Aerodynamics, and 
Controls) program (ref. 7). The present paper is the first which documents a major 
segment of the ISAC program. Future papers which document other segments may be 
found by employing ISAC as a key word in a library search. 

The equations of motion are outlined by employing a modal representation of the 
aircraft. Aerodynamic forces for either purely oscillatory or approximate fully 
unsteady motion can be considered. General linear control laws can be analyzed 



either by direct input of the elements of a transfer matrix or by construction uti- 
lizing selectable built-in filter types. 

Considerable effort has been directed toward the incorporation of features that 
allow the user to utilize his time efficiently. Interactive program execution is an 
option. Characteristic root variation with altitude, density, velocity, or control 
system gains can be performed in an automated manner with the results presented 
graphically . 

The STABCAR program and all user options are described in detail in this report. 
The functions of major program elements are specified, and key parameters are identi- 
fied that direct the input to, the computations in, and the results from each major 
subroutine. Flowcharts and a detailed description of input parameters are also 
provided. 

Results are presented and discussed which illustrate major program options. 
Characteristic root variations are shown with and without feedback control as a func- 
tion of parameters dependent on flight conditions. Root loci, as a function of vari- 
ation in control parameters, are also shown; comparisons of characteristi c roots 
obtained by using approximate fully unsteady and oscillatory aerodynamics are 
presented. 

The inputs for a series of sample cases are presented. The cases correspond to 
the graphical results that are shown. These sample cases provide explicit examples 
of the inputs required to exercise the major capabilities of STABCAR. 

STABCAR has been written specifically for efficient operation on CDC® CYBER 175 
equipment. Revisions are required to run this program on other equipment. 


OVERVIEW 

In STABCAR, stability characteristics are investigated. The dynamic elements 
that may be included are shown in figure 1. Mathematical descriptions of each of the 
dynamic elements are presented in the next section. The element representing the 
aircraft includes rigid, elastic, and control degrees of freedom as well as unsteady 
aerodynamic forces. The unsteady aerodynamic forces cause the equations of motion in 
the Laplace domain to be transcendental functions of the Laplace variable. Multiple- 
input/ mu It iple -output control laws can be studied; these can include dynamic compen- 
sation expressible in terms of a transfer matrix. 

STABCAR has been developed to solve for characteristic roots of the stability 
matrix. Typically, a search is made for roots as a function of the variation of a 
parameter of interest such as altitude, density, velocity, or a feedback gain. This 
procedure is illustrated in figure 2. Since the stability matrix is a transcendental 
function of the Laplace variable, the characteristic roots are determined by an 
iterative procedure. Consequently, initial estimates of the characteristic roots are 
required. Methods employed in STABCAR for determining initial estimates are presented 
in a later section. 


MATHEMATICAL REPRESENTATION AND SOLUTION TECHNIQUES 

In this section, the mathematical description of STABCAR is outlined. The equa- 
tions of motion are discussed as well as alternate options as to the form that can be 
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chosen for the aerodynamic forces. A capability for interpolation of beam- or plate- 
type mode shapes to obtain sensor inputs is described. The manner in which the con- 
trol system is included is described, and the techniques for characteristic root 
determination and major options for output of information are detailed. 


Equations of Motion 

The equations of motion consider only perturbations from a level equilibrium 
flight condition. The perturbed aircraft motion is represented in terms of a trun- 
cated set of rigid and elastic modes. The resulting equations, expressed in terms of 
generalized coordinates, masses, stiffnesses, dampings, and aerodynamic forces, are 
essentially the same as those developed in reference 8. The homogeneous part of 
these equations, expressed in the La place domain, are given in the following 
equation: 


2 

— 

( sb\ ( \ 

Ms + 

g. (s ) K. . 
i ,11 

- 

+ K + * F C ( T HS (s > H + T „E 1S) H e) 

— 


5 = { 0 } 


(1 ) 


where 

M ij 


K ij 


element of generalized mass matrix, // z. (x,y) a(x,y) z.(x, y) dS 

S 1 3 

element of generalized stiffness matrix; K(x, y;u,v) is influence 
coefficient giving force at (x,y) due to unit displacement at (u,v), 

jj z.(x,y) jj K(x,y;u,v) z.(u,v) dS dS 
c 1 o 3 \ Z 


sb 


-Q. . N ,— 
i j \ Ma U 


element of generalized aerodynamic force matrix with 

AP.(x, y;N ,p) being the pressure difference due to 
i Ma 

AP i (x,y;N Ma ,p) 

motion in jth mode, JJ z. (x.y) — - dS 

s 1 q 


g i (8) = ^ Tir g s, 


or 


If 


0 ) = 0 , 
n . 

i 


S . 0) 

i n. 
set to zero 


-? user chooses desired form. 


9 St structural damping coefficient associated with ith mode 

[ThsO] ng X n^ matrix of transfer functions relating actuator hinge 

moment outputs to sensor inputs 


t he ( s) 


n. x n diagonal matrix of transfer functions relating actuator hinge 

o 6 

moment outputs to actuator elongation 


n x matrix of modal coefficients converting hinge moment outputs to 

^generalized force 


5 

6 


vector, n x 1, of rigid and elastic generalized coordinates 
vector, n„ x 1, of control-surf ace generalized coordinates 


n c x i , of control-surf ace generalized coordinates 
6 
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z^(x,y) ith mode shape, z(x,y,t) = ^ z. (x,y) (t) 

i = 1 1 

H^j(x,y) linear or angular modal deflection contributing to input to sensor i 

resulting from unit deflection in jth generalized coordinate, 
r = H£ where H is n x n_ 

s 5 

H e n^ x matrix of modal coefficients relating actuator elongation to 

generalized coordinates 

An alternate form of the equations neglects aerodynamic hinge moments and hinge 
moments due to inertial coupling in which case the control surface deflections are 


6 = T (s ) H£ 
o 


The dimension of 6 is n^ * 1 and Tg is a matrix of transfer functions relating 
control deflections to sensor inputs* This expression replaces the control rows of 
equation (1); the resulting expression is 


M s^ + 


g i (s) hi 


+ + qQ cd N Ma'U 


sb 


— T ( S ) H 
o 


M r p 6 S + \Ma 9 U 


( 2 ) 


Unsteady aerodynamic forces *- In equations (1) and (2), the aerodynamic forces 
are expressed as functions of Mach number and -the complex variable p = sb/U. The 
programs currently available for production generation of aerodynamic forces can only 
compute these forces for p = \J~^1 ba)/U (e.g*, refs. 9 through 12). Two approaches 
are taken in STABCAR to circumvent this difficulty* 

The first approach, known as the p-k approximation (refs. 3 through 5) is valid 
for lowly damped modes. Two options are implemented for this approximation: 

(a) Form of reference 3 


Q(p) = Q(sb/U) » Q(p = 0 + \pt k) = Q(k) 


(b) British form (refs. 4 and 5) 

Q(p) = Q(sb/U) * Re(Q(k) ) 4- pQ* (k) 


Q* (k) is obtained by interpolation of the matrix of functions 


If k 1 = 0, 


r , /Q<V\ 

{o' (It 1 ) , Q’ (k 2 ), Q'(k n ^)} where O’ (k^ ) = Imf — j . 

Q* (k 1 ) = Q'(k 2 )* The approximation becomes exact at the flutter point. The second 
approach, known as the p-p method (ref. 3), requires the aerodynamic forces as a 
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function of p. In STABCAR, use is made of the concept of analytic continuation to 
develop an approximation to the aerodynamic forces for arbitrary motion in terms of 
knorn aerodynamic forces for oscillatory motion (p = ^ -1 k). The form assumed for 
the jth column of the aerodynamic force matrix is 


Aj 


(p) 


= A. 


+ A 1 p + A 2 p 


+ E 


A U+ 2) , P 
1 


p + b 


Aj 


A=1 


(3) 


This expression has the same form as that employed in references 1 and 2. Determina- 
tion of the polynomial coefficients is initiated by computation of k) for a 

number of values of reduced frequency. Ihis computation must be done^by an external 
program such as those described in references 9 through 1 2. The user specifies the 
vector bj and may impose linear constraints upon the coefficients in equation (3). 
For example, the Aq columns associated with, rigid-body linear displacements should 
be zero. The coefficients of each element of ^ (p) are determined, subject to any 
imposed constraints, such that the errors between the curve fit (eq. (3)) and the 
tabular data are minimized in a least-squares sense. The leas t-squares -fit procedure 
and allowable constraints are described in appendix A. 


References 13 and 14 present a theoretical basis for the curve-fit procedure and 
describe the asymptotic form Q. (p) should take. The asymptotic limits are not 
enforced in STABCAR because constraining the curves to fit high frequency limits 
tends to degrade the fit in the limited frequency band where tabular data are avail- 
able. Reference 1 5 describes the application of a constrained optimization procedure 
to optimally choose bj in equation (3). 

References 16 and 17 describe a program now under development which will allow 
computation of aerodynamic forces for transient motion where p has a nonzero real 
part. STABCAR can be readily modified to utilize such data. The required modifica- 
tion is to provide for interpolation of tabular aerodynamic force arrays as a func- 
tion of two independent variables (the real and imaginary parts of p) rather than as 
a function of only the imaginary part of p. 

Sensor input definition .- Sensor inputs in STABCAR are assumed to be linear 
combinations of the generalized coordinates. For linear displacements. 


r D (x,y) = ^ z_.(x, y) £. 

j = 1 


For angular displacements. 


r Q (x, y) 


-E 

j=i 


5z . (x, y) 


dx 
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y 

These expressions would be multiplied by s for a velocity type sensor and by s 
for an accelerometer. The multiplication by the appropriate power of s is included 
within the computation of the effects of sensor dynamics. (See eq. (8).) Thus, in 
computing sensor inputs, the assumption is made that inertial axes are employed. The 
user may override this assumption by appropriately modifying the rigid-body portions 
of the mass, damping, and stiffness matrices. 

The mode shapes or their derivatives at the sensor locations must be specified 
before the sensor inputs can be determined. Two options are included for defining 
the modal data for input to sensors. The first option is by direct input of the 
modal data required at the sensor locations. The second option is to input all the 
modal data (mode shapes and node-point locations) and interpolate for the modal 
amplitudes at the sensor locations. This option is attractive when performing 
sensor-location studies. 

The remainder of this section describes the interpolation procedure. Vibration 
analyses, which use modelings ranging from simple beams to finite-element representa- 
tions, calculate modal data at a number of grid points. In order to accommodate out- 
put from a variety of vibration analyses, it is necessary, in STABCAR, to arrange the 
modal data into subset divisions called structural sections. For aircraft configura- 
tions, the structural sections follow naturally as the wing, fuselage, horizontal 
tail, and so forth. 

There are two types of interpolation methods in STABCAR corresponding to the two 
types of structural sections available in the program. The first is a displacement- 
twist interpolation method which is most commonly used when the mode shapes are gen- 
erated from an elastic-axis beam model. The second is a surface-spline interpolation 
method (ref. 18) which is used when the mode shapes are generated from a finite- 
element "plate-type " model. The entire modal representation may consist of a combi- 
nation of both types of structural sections. 

For a zero-thickness idealization of a surface and a beam structural representa- 
tion, the displacement in the ith mode normal to the surface can be written as 


z . (x, y) = z . (x* , y 1 ) = z . (y ' ) 

li l 


x'R , (y ' ) 

y’ 


(4) 


where 


(x, y,z) right-hand coordinate system with x positive aft and z positive up 

(x'jy'jZ 1 ) right-hand coordinate system with z 1 parallel to z and x* 
and y* swept from (x, y) by angle A 

z^(x,y) deflection in ith mode at (x, y) normal to surface 

z^(y') vertical bending in ith mode at (x' = 0,y') 

R ( y 1 ) rotation about y* axis in ith mode at (x 1 = 0,y') 

\T * 
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Figure 3 shows the geometrical relationships for a flat surface parallel to the xy 
plane. The origin of the swept, translated axis system is at (x 0 ,y 0 ); A is the 
sweep angle (positive clockwise); and positive z is out of the page. The transfor- 
mation from (x, y) to (x^y 1 ) is given by 


fl 

u 

cos A 

-sin A 

( x - x o| 

(y'J 

i 

sin A 

cos Aj 

- y <j 


For angular displacements, the partial derivative of equation (4) is required 
and is given as follows: 


9z . 
x_ 

5x 


dz . (y 1 ) 
i 


ay' 


- x 1 


dR _ (y f ) 

y l 


dy’ 


sin A - R , cos A 

y i 


or 


dz . 

i_ 

9x 


R , <y f ) 

x! 

i 


- v' 


dR y , (y 1 ) 


dy ' 


sin A - R . (y * ) cos A 

y i 


( 6 ) 


where R f is rotation about x 1 (positive counterclockwise), 
i 

For !! plate !! type modes the modal data from a structural analysis are defined at 
points distributed over the structural section. Thus, there are two independent 
variables (x,y) and the i nterpolation is done by using the surface spline methods 
described in reference 18. Slope information (e.g., 5z,(x,y)/5x) is generated by 

interpolating for the derivative of the surface spline. 

Rigid-body mode shapes may be defined by input by the user. In this case, the 
rigid-modal data must be treated as if they are elastic modes (i.e., n r = 0 where 
n r is the number of rigid modes). Alternatively, if the symmetric rigid-body mode 
shapes of vertical translation and pitch about the aircraft center of mass are 
desired, the user may set n r = 2. Then the following built-in expressions are 
employed to compute the rigid-body contributions to sensor inputs: 

linear displacements: 


z (x,y) = 1; z (x,y) = -(x - x )0 
i 2 eg n 


7 


Angular displacements: 


dz 

dx 


<x,y) 


= 0 ; 


dz 


dx 


-0 

n 


Control-System Representation 

"The control system is divided into three parts: sensor, compensation (filter- 

ing), and actuator dynamics. The input-output relationships between these dynamic 
elements and their coupling with the aircraft dynamics are indicated schematically in 
figure 1 and explicitly in equations (1) and (2). In a subsequent discussion of 
these dynamic elements, the term "block" is used interchangeably with "transfer 
matrix. " 


Sensor dynami cs options .- Input to the sensor dynamics can be either linear or 
angular position information expressed in terms of the generalized coordinates as 


r - H£ 


(7) 


The dimension of r is n x 1 . Within the sensor dynamics block each r. can be 
transformed to a position, velocity, or acceleration input by multiplication by the 
appropriate power of the Laplace variable s. A given sensor output will be of the 
f orm 


r * 

D 


nSn .-1 

a^ + a_ s + • . . + a s ' 

Sn . . Sn . Sn . „ 

1 j 2j (nSn . ) j J . 

J J 

=> s r . 

nSd . -1 l 

a e . + a s + ... + a s 3 

ID 2d (nSd )] 


T sj (s) 


r . 
J 


( 8 ) 


where Jj =0, 1, or 2 corresponds to a displacemen t , rate, or acceleration sensor, 
respectively. Default values are built into the program so that only Jj must be 
input if sensor dynamics are neglected. 

Compensation options .- The transfer matrix T^, of dimension n^ x n g , specify- 
ing corrpensation of the signals from the sensors prior to sending them to the actua- 
tors can be full with distinct dynamics in each element. The transfer matrix rela- 
tionship between compensator outputs and compensator inputs (sensor outputs) is 


6 = T (s) r 1 

Ij Li 


( 9 ) 


where 
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Here (G)^. is a gain, each [t l ( s ) j may contain dynamic compensation elements 
but has a default value of 1. Alternate values for [T L (s)]j_j may be defined either 
by direct input of polynomial coefficients or by utilization of the following 
built-in filter types: 

1. Notch filter: 


1 4 - 2s 




where u) and £ are the natural frequency and damping ratio of the notch filter 


2. Integral filter: 


V s 

3. Proportional plus derivative filter: 

K 1 + K 2 s 

4. Lead-lag (lag-lead) filter: 

1 + S 

1 + x s 
2 


5. Polynomial filter: 


+a„ s + . . . + a^ s 

Fn . Fn ^ . Fn . ^ x . 

1 1 2i (nFn . )i 


nFn . -1 

l 


a P , + a„_ s + ... + a s 

1 1 2i (nFd . ) i 


nFd , -1 
x 


(1 0 ) 


Any combination of these filter types may be utilized including repetitive use of the 
same filter type with the same or distinct constants. Additional filter options are 
as follows: 

6. Scheduling: 

The capability for providing scheduled compensation elements T LS (s ' N Ma' q) ij 
is also included. For each control-sensor pair, it may be desirable to 
allow the control law to vary as a function of Mach number and dynamic pres- 
sure. This variation is accomplished by providing a call to a subroutine 
which has available for use control-sensor pair identification, the value of 
the Laplace variable s, Mach number, and dynamic pressure. The subroutine, 
which defines the desired scheduling, must be supplied by the user. This 
procedure is illustrated in appendix B for a specific scheduling law. 
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7. Phase error: 


An additional feature included is the following option of introducing a 
phase error , into the compensation for one control-sensor pair; 



This option allows one to study sensitivity to phase errors. 

Actuator dynamics options .- The transfer-function matrix relating inputs to and 
outputs from each actuator is diagonal. Its default value is the identity matrix. 

If actuator dynamics are to be considered, the following numerator and denominator 
polynomial coefficients are input: 


(6 ). = 

A l 


a +a s + . . . + a £ 

An . An ^ . An . . . 

1i 2i (nAn . )i 

l 


+ a s + . . . + a s 

Ad . Ad . Ad, _ 

1 1 2i (nAd . )i 

l 


nAn . -1 

l 


'nAd. -1 (6 lA 
1 


(1 2 ) 


or 


ft)’ [\]' 


where the dimension of 6, is n c x i. 

A 6 

Combination of sensor, compensator, and actuator matrices .- The relationship 

between 6 and r is 
A 

If aerodynamic and mass coupling hinge moments are considered, the resulting transfer 
matrix corresponds to T HS in equation (1); otherwise it corresponds to in 

equation (2). The polynomial d(s) in equation (13) is the common denominator of 
T^(s) (see eq. (9)). This factor can be cleared from the denominator in equa- 
tion (1) or (2). Appendix C describes this option. 


Characteristic Root Determination 

The primary objective of STABCAR is to determine characteristic roots of an 
aircraft represented mathematically by the transcendental matrix equation (eq . (1) 

or (2)). Roots can be found as a function of variations in velocity, density, alti- 
tude, or control parameters. 


Determinant iteration .- Characteristic root determination is accomplished mode 
by mode with determinant iteration given initial estimates for the roots. The roots 
to be traced are specified by input. Methods for obtaining the initial estimates are 
described after the following discussion of the procedure for determining a charac- 
teristic root. 

An attempt is made to find the ith characteristic root in a manner which avoids 
convergence upon a previously determined root, j < i. This is accomplished by 
determining the zeros of 


D' (s) = D(s) (i = 1) (1 4a) 

or 

D' (s) = , - P(s) (i > D (14b) 

1-1 

n (s - x. ) 
j-i : 

where 


D(s) value of complex determinant of stability matrix in equation (1) or (2) 

X. previously determined characteristic root 


Let + 1 denote the desired n + 1 estimate of the zero of equations (1 4) 

corresponding to the ith characteristic root and denote a previously determined 


estimate; therefore, u : 

i 


or d? is the value of D* (s) at 


n+1 


L n+1 


or 


respectively. There exists an a > 0 such that, for 

n 


-n + 1 


= ss 


a D? 
n i 


dD! 

i 


Ids 


(1 5) 


d; 

< 

D! 

i 


l 

n+1 


n 


In the search for Xj_, a is constrained to be <0.1 in order to reduce the 
likelihood of jumping to a point sufficiently far away from X. as to cause converg- 
ence upon another root. The general procedure is to solve equation (1 5) repetitively 

until si differs from s-; by less than a tolerence e specified by the user. 

n + 1 -Ln 


1 1 


That is, if 


Si n+1 



< e 




(1 6a) 


or 


Si n + 1 




< e 



(1 6b) 


convergence has occurred and Ka - s\ . • 

n+ i 

Initial estimates for roots .- The algorithm for determining characteristic roots 
described in the previous section requires both function (eqs. (14)) and slope 
(eq. (1 5) ) information. Thus, two initial estimates are required for each root that 
is to be traced. Three options are available for specifying the initial estimates: 

1. Direct user input: 

2. Computations based upon system natural frequencies of vibration in vacuum: 


s i, = 0 + (-i <% 

s i2 = (d! + \TT d 2 )w ni 






J 


(1 7) 


where d^ and d 2 are input parameters having default values of 0.01 and 
1.0, respectively. This approximation works well for the elastic degrees of 
freedom at sufficiently low dynamic pressures. 

3. Initial estimates by matrix iteration: In this approach the matrix equa- 

tion (1) or (2) is expressed in the form 


s 2 l = A(s) l 


(1 8a) 


or 


2 

s 



= A(s) 



(1 8b) 
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Appendix C shows how A(s) is obtained, and appendix D describes the matrix itera- 
tion method* 

Prediction of and 2 at a new va lue of independent variable.- Subse- 

quent predictions for s-j^ and s-^ at each new value of the independent variable 

are obtained by using an extrapolation of either a linear fit to the two preceding 
converged roots (after second iteration), a quadratic fit to the preceding three 
converged roots (after third iteration), or a cubic fit to the preceding four con- 
verged roots (after the fourth iteration). The type of extrapolation to be performed 
(linear, quadratic, or cubic), after the requisite number of iterations, is specified 
by input. 

Logic during convergence difficulties .- Failure to converge upon a characteris- 
tic root occurs periodically when the root is changing rapidly as a function of the 
parameter being varied. When such a failure occurs, the matrix iteration approach is 
employed in an attempt to recover. The matrix iteration approach is started with the 
last successful estimate for the root. If this attempt is successful, normal deter- 
minant iteration is then resumed. If the attempted recovery by using matrix itera- 
tion fails, the following scheme is utilized: 

1. Let \A denote the value of the parameter being varied (independent 

variable) at the point the root was last successfully obtained. Let V^ + 1 

denote the value of the independent variable at which convergence failed. 

2. Divide the interval (V^ , ) into a user specified number of increments 

n 

AV 

c 

3. Solve for the root by using determinant iteration at each intermediate point 

and the end point + 1 . 

If convergence is achieved at \L +1 , normal execution is resumed. If not, that par- 
ticular root is deleted from the set being sought and normal execution is resumed for 
the remaining roots. 

Precise determination of point at which change in stability occurs .- When the 
first change in stability is noted, the value of the parameter being varied which 
corresponds to zero for the real part of the root is automatically determined. The 
determination is accomplished by 

1. Dividing the interval between where a change in stability is noted, 

and V ± into a user specified number of increments n^ v 

v s 

2. Evaluating the characteristic root at all intermediate points 

3. Performing an interpolation with these data to find the precise value of V 

at which the real part of the characteristic root is zero 

Matched-point computations. - Two options for automatic matched-point computa- 
tions are provided. Each seeks matched points at a particular input Mach number. A 
matched point is one for which the density, Mach number, and velocity are consistent 
with actual physical conditions (e.g., atmospheric conditions or tunnel operating 
conditions). Characteristic roots are computed and their corresponding root loci or 
damping ratios and natural frequencies can be plotted as a function of the parameter 
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being varied. The point at which the real part of a root becomes zero denotes a 
change in system stability. For an elastic mode, this corresponds to a matched flut- 
ter point. The two automated types of matched-point computations are 

1. Wind-tunnel matched point (density variation at fixed Mach number): 

(a) Characteristic roots are examined to determine stability characteristics 

with velocity held fixed. 

(b) Outputs for this type run include damping ratios and natural frequencies 

or root loci as a function of either density or dynamic pressure. 

Density variation analyses correspond to the method of operation of some wind 
tunnels and are useful for correlating analytical and experimental results 
for scaled flutter models. 

2. Atmospheric matched point (altitude variation at fixed Mach number): 

(a) Beginning at a high -altitude (low-q) condition, characteristic roots are 

determined at specific altitudes with the data of reference 19 to 
determine the corresponding density and speed of sound. This type run 
yields a matched flutter point, if one exists, for the specified Mach 
number. 

(b) Oatputs for this type of run include damping ratios and natural 

frequencies or root loci as a function of either altitude or dynamic 
pressure. 


PROGRAM DESCRIPTION 

In this section, program features are outlined and the primary functions of 
major program modules are identified. 


Features for Enhancement of Resource and User Efficiency 

Considerable effort was expended to develop a code that makes efficient use of 
resources and allows the user to efficiently study the stability characteristics of 
flexible aircraft. The following features are a result of the effort. 

Overlay structure .- STABCAR is an overlay program with two levels of overlays. 
The overlay structure is indicated in figure 4. Overlaying allowed modularization of 
the code which facilitates modifications affecting only one overlay. It also reduced 
the amount of code required in core at one time. 

Variable dimensions .- All large arrays are variably dimensioned, with the dimen- 
sion specified by input. This enables the user to either study a new larqer problem 
or make efficient use of central memory when a smaller problem is to be investigated 
without making program modifications. 

Dynamic storage allocation .- Within each overlay, the required computer memory 
is computed and updated at critical points during job execution. This dynamic allo- 
cation of storage results in more efficient utilization of central memory. Further 
description of this technique is given in reference 20. 



Free-field reads ,- Free-field reads are employed to input most of the array data 
(e.g. , mode shapes, generalized masses, stiffnesses and dampings, and unsteady aero- 
dynamic forces). 

Model selection .- A subset of the model that has been input can be selected for 
study. This option allows one to consider specified degrees of freedom while ignor- 
ing others that are represented in a given set of data. The appropriate generalized 
masses, stiffnesses, aerodynamic forces, and so forth are automatically extracted. 

Automated parameter variations .- Stabi lity characteristics are typically studied 
as a function of variation of dynamic pressure or control-law parameters. Automatic 
variation of dynamic pressure can be specified by the user in three ways: (1 ) alti- 

tude variation, (2) density variation, or (3) velocity variation. Gain variations 
can also be performed in an automated fashion. The effects of variations of other 
parameters - for example, sensor locations - upon system stability can be examined by 
running several cases, each having a different value of the chosen parameter. 

Graphical display of output. - Stability characteristics can be displayed graphi- 
cally in two forms. The' damping ratio and natural frequency corresponding to each 
characteristic root can be plotted, and characteristic root loci can be plotted as a 
function of any parameter which can be varied in an automated manner. In addition, 
plots can be made which show the variation of the unsteady aerodynamic force with 
reduced frequency. The aerodynamic-force plots can exhibit the data points input 
into STABCAR, curves resulting from interpolation, and curves resulting from a 
p-plane approximation to the input data. The Langley Graphics package employed to 
obtain graphics output is not included in STABCAR. A description is given in appen- 
dix E of the function of the routines in this package to facilitate substitutions of 
equivalent software. 

Interactive execution .- Interactive execution of STABCAR is an option. Some 
advantages of interactive execution are the more rapid receipt of output, the immedi- 
ate visibility of input errors, and the capability of early termination. 

Restart capability. - Data that define the mathematical model being employed in 
the current case are stored on TAP El, a random access file. The data include 
selected modal information and corresponding aerodynamic forces including a p-plane 
fit if it has been generated. In addition, if the plot control parameter is nonzero, 
characteristic roots are stored as a function of the parameter being varied along 
with the user-supplied case title for identification. Consequently, if TAPE1 is 
saved, it may be used in a subsequent case if the model selected is a subset of the 
previous case. Thus, for example, the p-plane coefficients need not be recomputed. 
The new case might be to extend the range of variation of a parameter, to delete 
redundant or erroneous characteristic root data, or to combine sets of characteristic 
root data into one composite plot. Examples are given in the sample cases which 
illustrate several of these possibilities. 


Functional Description of Overlays 

The program is structured with two levels of overlays. The first level (overlay 
MAIN) is an executive. A brief summary of the major functions of each overlay will 
be presented here. A more detailed description including the function of major sub- 
routines and definition of key parameters which direct the computations is given in 
appendix F. 
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MAIN: Ibis executive overlay directs the computations by calling other overlays 

based upon input parameters. Namelist data are input which define the type of stud- 
ies to be performed. 

INTERP: INTERP performs geome try -re la ted computations needed for a surface 

spline fit if mode-shape data are for a plate. The data required in these computa- 
tions are input here. The results of the computations are subsequently used in 
SENSOR to generate surface spline coefficients that fit the plate modal data. 

SENSOR: This overlay determines surface spline coefficients for fitting plate- 
type modes. Beam-type modal data are also input here. Interpolation of beam and 
plate modal data is performed to obtain modal coefficients defining either deflec- 
tions or slopes that would be measured by sensors at specified locations. These 
modal coefficients are stored for subsequent use and define the matrix H in equa- 
tion ( 1 ) or ( 2) . 

MATINPT: Large data arrays such as the aerodynamic data, generalized masses, 

generalized stiffnesses, and control-system dynamics definitions are input. In addi- 
tion, the selection feature is contained, which enables one to choose a mathematical 
model that is a subset of the model defined by a given set of data. 

PPLANE : This overlay is called if a p-plane approximation to the aerodynamic 

forces is desired. The best coefficients are determined, in a least-squares sense, 
to employ in equation (3), given a set of aerodyamic force data defined over a range 
of reduced frequencies. These coefficients are stored for subsequent use in overlays 
PKFLUT and AEROPLT. 

CONTROL: The matrix of transfer functions T T * T is constructed. (See 
eqs. (9) and (13).) The ( )— element relates tfte^output of actuator i to the 
input of sensor j excluding (G)^ and any scheduling or phase error. The poly- 
nomial coefficients for sensor, compensation, and actuator dynamics are constructed 
separately and then combined. Each element of the resulting matrix can then be mul- 
tiplied by a distinct feedback gain and scheduling component in overlay PKFLUT. 

PKFLUT: The matrix in equation (1) or (2) is formed. Characteris tic roots are 

determined as a function of variation in altitude, density, velocity, or feedback 
gains and the resulting output is stored for possible subsequent plotting. By oper- 
ating in a multicase mode, the effect of variation in other parameters such as actua- 
tor or compensator dynamics or sensor location can also be studied. 

AEROPLT: Plots are constructed which show hew the oscillatory aerodynamic 

forces vary with reduced frequency and depict how the p-plane approximation fits the 
data. 


PKPLOT : Plots showing root loci resulting from variation in velocity, density, 

altitude, dynamic pressure, or feedback gains can be generated. Alternatively, damp- 
ing ratios and natural frequencies can be plotted for specified modes as a function 
of any of these parameters. Namelist inputs are accepted which define the types of 
plots desired; scales; plot titles and which plot data files, among those available, 
should be combined to make a composite plot. 



STABCAR Flowchart 


A flowchart is presented in figure 5, which describes the computational flow in 
STABCAR. This flowchart, although relatively detailed, is not intended to describe 
either the effects of all input options or the function of all program subroutines. 

INPUT REQUIREMENTS 
Description of Input Files 

In this section STABCAR input requirements are defined. The documentation of 
the input requirements is presented in tables I through V. Therein, the types of 
data are outlined; a detailed description of the inputs is given which defines the 
individual input variables along with case -dependent conditions that determine 
whether certain of the data are required? and the input files that contain the data 
are identified. The correspondence between the FORTRAN names and the algebraic names 
of the variables is also indicated in these tables. 

TAPE2 contains data which define problem size, user options, starting estimates, 
control system definition, and so forth. These data, primarily namelist, are out- 
lined in table I. Table II contains a detailed description of the TAPE2 inputs. If 
ICASE is positive, the run is to be an interactive one where changes in data are made 
on a remote terminal during program execution. In this case, all such changes in 
input are made interactively from file INPUT during program execution. 

TAPE5 is a file containing the bulk of the array data required for program exe- 
cution. All such data are read free field. The data contained on TAPE5 and the 
array sizes are shown in table III. 

TAPE1 0 contains mode-shape data to be employed in determining sensor deflections 
for computing feedback signals as indicated in table IV. 

TAP El is a random-access binary file created by the program. Data on TAP El have 
been either input or generated in a particular overlay of STABCAR. These data may be 
required in other sections of the program during the current case or may be needed in 
a subsequent case. For example, the user may have selected a subset of the mathe- 
matical model. This file could be saved by the user for future runs. The option of 
using the TAPE1 data as input is accomplished by setting ICASE = ±1 and the option 
of selecting a subset of the model data stored on TAPE 1 is exercised by setting 
I SELECT = -1 as illustrated in case 6 in appendix G. The contents of TAPE1 are 
indicated in table IV. 


Input Aids 

Many cases are run which employ a previously generated TAP El for the bulk of the 
input. The user may be in doubt as to precisely what values are stored on TAP El . 
TAPE2 variables and selected sensor data stored on TAPE1 may be determined by execut- 
ing STABCAR with IFLUT = 0, ICASE = ±1, ICONSYS = 1, and TAP El local. If IAEROWR 
is set equal to 1, the aerodynamic force tabular arrays will also be output. If 
IMODEWR is set equal to 1, generalized masses, structural dampings, and either 
natural frequencies or the generalized stiffness matrix will be output. 
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A local scratch file called STABIN is produced by STABCAR which reflects TAPE2 
or INPUT file data generated during a series of batch or interactive STABCAR runs. 
When ICASE = ±2 (i.e., an ini tialization run), STABIN contains the values of all 

variables including default values, except for array variables where only the first 
value and any nonzero values are exhibited. When ICASE = +1, STABIN contains the 
values actually input as well as any default values pertinent to the analysis being 
attempted. STABIN is not rewound by STABCAR; therefore, it can provide the user with 
a record of what was done during a multicase session. 

When IPKLT = +1, the user has a number of options pertaining to any sets of 

plot data stored on TAPE1 . If the user sets IDELPLT = +1 , he may delete sets of bad 

plot data by specifying which of the stored plot sets he wishes to keep. The user 
can specify which sets he wishes to keep in any order. They will then be resequenced 

in the order selected. If IDELPLT = 0, the user may combine a series of plot sets 

to form a composite plot. (See appendix G, table GVII, for example). Here too, the 
plot sets may be chosen in any convenient order. The choice of order is generally 
made so as to maximize ISAME (table II), which minimizes the number of times namelist 
$PL0TP has to be input (table GVII). 


PROGRAM CAPABILITY DEMONSTRATION 

In this section some of the capabilities of STABCAR are demonstrated. A mathe- 
matical model of an aircraft is defined and STABCAR is employed to examine its sta- 
bility characteristics as a function of several parameters. The inputs required to 
generate each sample case are shown and a graphical display of the results from each 
case is presented. 


Aircraft Mathematical Model 

Results obtained by using STABCAR are shown for a mathematical model of an early 
version of the DAST ARW-2 (Drones for Aerodynamic and Structural Testing, Aeroelastic 
Research Wing Number 2). The DAST ARW-2 has a high -aspect-ratio (10.3) supercritical 
wing having 25° sweep at midchord mounted on a Firebee drone fuselage. The struc- 
tural design of the wing was performed by the Boeing Military Airplane Company 
(refs. 21 and 22) under the assumption that an active flutter suppression system 
would be operational. The wing will flutter within the flight envelope with the 
flutter suppression system turned off. The structural model was supplied to NASA by 
the Boeing Military Airplane Company. The structural model was input into the 
NASTRAN program (ref. 6) to obtain the modal characteristics of the model. Two 
rigid-body modes (plunge and pitch), 10 symmetric elastic modes, and 1 control mode 
were retained for computation of generalized aerodynamic forces. The resulting modal 
information together with geometric and paneling data were input into the ISAC 
(Interaction of Structures, Aerodynamics, and Controls) program (ref. 7) where 
unsteady aerodynamic forces were generated by using a doublet lattice method. The 
aerodynamic paneling employed, as well as the location of the control surface and 
sensors to be employed for a flutter suppression system, is indicated schematically 
in figure 6. The flutter suppression control law that is used is shown in figure 7. 
In this control law, the difference in signals from two accelerometers located at 
(x s ,y s )i = (7.1 44,2.083) m and ( X g /V s ^2 = ( 7 • 3 36, 2 . 1 34 ) m is fed back to an actua- 
tor driving an outboard aileron. Note that there is a scheduled component in the 
filter which varies with Mach number and dynamic pressure to provide appropriate 
phasing between sensor output and actuator input at varying flight conditions. 
Appendix B shows how the scheduled component was incorporated into STABCAR. 



Graphical Output From Sample Cases 


The amount of information to be output by STABCAR is selectable by the user. 

The basic file for printed output is TAPE6 which can be sent to the printer for 
either batch or interactive runs. This lengthy output is well identified by 
Hollerith descriptors and is not shown here. When interactive runs are being made, 
the user will typically minimize the information displayed to the screen. The vari- 
ables employed to designate the output desired have been defined in table II in name- 
list $INPUT. In the results which follow, only graphical outputs to the screen are 
presented. The inputs required to generate each figure are given in appendix G. 

Case 1 - Aerodynamic forces and p-plane approximation .- In case 1 (aerodynamic 
forces and p-plane approximation), a p-plane approximation was made to the aerody- 
namic forces. Figure 8 shews the selected aerodynamic force data input into STABCAR, 
the curves resulting from a quadratic interpolation for points between the data, and 
a p-plane fit. The elements shown are those for the three elastic modes that are 
most critical in the flutter mechanism and those for the corresponding control sur- 
face effectiveness elements. The inputs required to generate these figures are pre- 
sented in appendix G. 

Case 2 - Open-loop characteristic roots as function of altitude .- Figure 9(a) 
shows the root loci of the aircraft without controls as the altitude is varied. 

Case 2 (open-loop characteristic roots as function of altitude) is for a Mach number 
of 0.86; each point is matched by solving for the atmospheric density at the altitude 
as well as the speed at that altitude corresponding to the 0.86 Mach number. The 
open-loop flutter occurs at h = 6.705 km which corresponds to a velocity of 
269.6 m/sec and a density of 0.6101 kg/m^. Note that the roots corresponding to 
modes 3, 5, and 8 are varying with altitude, whereas the other roots considered are 
relatively unaffected. It will be shown, in a subsequent plot, that essentially the 
same flutter point is predicted with a mathematical model that retains only the 3d, 
5th, and 8th modes in the analysis. Figure 9(b) presents essentially the same infor- 
mation as figure 9(a) except damping ratios and natural frequencies are plotted as 
explicit functions of the altitude. 

Case 3 - Open-loop characteristic roots as function of density .- Results 

obtained by running case 3 are shown in figure 10. Figures 10(a) and (b) show root 

loci and damping ratio and natural frequency variations, respectively, as a function 

of density for N w = 0.86 and U = 269.6 m/sec for no controls. The flutter point 
Ma o 

occurred at p = 0.6101 kg/m which is in perfect agreement with the altitude varia- 
tion case. 

Case 4 - Open-loop characteristic roots as function of dynamic pressure .- Fig- 
ures 11(a) and (b) for case 4 (open-loop characteristic roots as function of dynamic 
pressure) correspond precisely to figures 10(a) and (b). They are density variation 
runs but the independent variable has been changed to dynamic pressure to show the 
characteristic roots as a function of q. Plots of root variations as a function of 
dynamic pressure can also be made in altitude and velocity variation runs. 

Case 5 - Characteristic roots as function of velocity (truncated model) .- Case 5 
shows results from using a truncated model. Figure 12(a) is an open-loop plot of 
damping ratios and natural frequencies versus velocity. The aerodynamic forces for 
= 0.86 and the matched-point density found in figures 9 and 10 are employed. 
Thus, only the point at which U = 269.6 is a matched point. In this run, the 
select feature has been employed to consider only the effects of the modes that are 
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most important to the onset of flutter. These are modes 3, 5, and 8 which are the 

1st, 3d, and 6th elastic modes. The predicted flutter speed of 271.2 m/sec is within 
1 percent of that found earlier with two rigid-body and six elastic modes. In fig- 
ure 12(b) damping ratios and natural frequencies are plotted as a function of veloc- 
ity under the same conditions and assumptions as for figure 12(a) except that the 
control loop is closed. The control-sytem representation is that shown in figure 7. 
The aircraft does not flutter for the range of velocities shown. 

Case 6 - Elastic mode gain root loci. - Figure 13 shows gain root loci (case 6) 
for elastic mode roots at = 0.86 and h = 4.5 72 km which correspond to U Q . 

For this condition the open-loop ARW-2 configuration is unstable as indicated by the 
root in the right half -plane in figure 13. The root loci are shown for both a p-k 
(solid lines) and a p-p (dashed lines) analysis. The unstable mode is stabilized at 
a gain of approximately 0.7. The slight difference between the p-k and p-p analyses 
for this root at the ju)-axis crossing is an indication that the least-squares fit of 
the aerodynamic forces is not precise at the flutter frequency. Note that the mode 
number identification does not correspond to the numbering of the earlier figures. 

The matrix iteration technique was employed here to obtain initial estimates for the 
open-loop roots. The roots are numbered in ascending order of magnitude of the open- 
loop roots. Since the two rigid-body roots are not shown, the lowest number is 3. 

Case 7 - Effect of phase errors upon elastic mode gain root loci. - Figure 1 4 
illustrates the use of STABCAR to determine whether specified phase and gain margins 
are met by a candidate control law. This figure presents the gain root loci of the 
elastic mode roots (case 7) at the same flight condition as those for figure 1 3 for 
the p-p analysis only. Additional root loci are shown where the control signal has 
been modified by phase errors of -45° and 45°. The figure indicates that the candi- 
date control law has less than a -45° phase margin at the nominal gain and does not 
meet a -6 dB gain margin at nominal phase. 

Case 8 - Gain root loci of control system poles .- Figure 15 shows the gain root 
loci for the control system poles (case 8) at N Ma = 0.86 and h = 6.705 km, which 
is the open-loop flutter point. The dynamic pressure at this condition is 25 percent 
less than that at U D (the flight condition of figs. 13 and 14). For this dynamic 

pressure, root loci are shown for both p-k (solid lines) and p-p (dashed lines) 

analyses. Loci on the negative real axis were determined for the p-k analysis only 
because of the complications occurring in the p-p analysis due to isolated singular 
points and poles introduced by the p-plane approximation to the aerodynamic forces. 
Note in particular the branch where two real poles merge to form a complex pair. 

This type of locus can best be handled in an interactive run by trying different 

initial estimates having nonzero imaginary parts when the branch point is approached. 
Another curve of interest is the locus of the filter root which has an open loop pole 
at (-50,1 96.5) rad/sec. Its locus is similar to that of the 3d elastic mode in fig- 
ure 13 which emphasizes the variation in root loci patterns that can occur as a 
result of different zero and open-loop pole locations between two flight conditions. 

Case 9 - Short-period altitude root locus .- A short-period root locus (case 9) 
with altitude is presented in figure 16 for p-k and p-p analyses. The rigid-body 
degrees of freedom considered in this study are plunge and pitch about the center of 
mass with respect to an inertial set of axes. Drag contributions to the perturbation 
equations are neglected. The resulting equations for the rigid-body roots correspond 
approximate ly to those for the short-period approximation with additional effects of 
unsteady aerodynamic force contributions due to both rigid and elastic motion. 
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CONCLUDING REMARKS 


STABCAR is a computer program designed for the study of the stability charac- 
teristics of flexible, actively controlled aircraft. Analyses can be performed that 
include either oscillatory or approximate fully unsteady aerodynamic forces. Deter- 
minant iteration is employed for characteristic root determination. Consequently, 
sufficiently accurate initial estimates and sufficiently small steps in the quanti- 
ties being varied are required to converge upon and follow the characteristic roots. 

STABCAR can be executed interactively with a remote terminal which allows the 
user to rapidly assess the results he is obtaining and conveniently change the input 
data. Changes such as deleting modes from the model, altering transfer functions or 
sensor types, or executing different types of analyses may be performed quickly and 
easi ly. 

In an effort to aid potential users, a description of the mathematical tech- 
niques employed has been outlined, a flowchart has been provided, key variables have 
been defined, and a description of the input requirements has been given. In addi- 
tion, a series of sample cases have been presented which provide explicit examples of 
the inputs required to utilize the major program capabilities. 

STABCAR has been written specifically for efficient operation on Control Data 
equipment. Revisions are required to run this program on other equipment. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
March 7, 1984 
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p-PLANE APPROXIMATION TO THE UNSTEADY AERODYNAMIC FORCES 
Least-Squares Solution 

It is desired to approximate the jth column of aerodynamic forces with the fol- 
lowing function: 


Qj (P) = A 0. 


A 1jP + A 2 jP 



A U+2) jP 
P + 


(A1 ) 


The form of this expression is the same as that found in references 1 and 2. How- 
ever, constraints would have to be imposed upon the coefficients to have full equiv- 
alence between the fit described here and that of references 1 and 2. 

Assume tabular data for oscillatory motion have been computed: 




Q. OO* 

D 2 



(A2) 


The objective is to determine the coefficients in equation (A1 ) such that 

qlj (p = k) best fits the tabular data in a least-squares sense subject to a set 

of linear equality constraints which may be imposed upon the coefficients. 


Consider a particular element 


2 i j 


of 




Let 



A (, 


Aj 



(A3) 


denote the coefficients to be determined. Suppose that there are linear equality 
constraints which the coefficients must satisfy denoted by 


T 

c;a. . 
3 ID 


C 


0 . . 

ID 


The error in the fit of equation (A1 ) to the data at k - k^ is 



r 
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j 
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L il 
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Re f Q. . 
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13 

n 

> 


il 

n J 


(A4) 


(A5) 
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or 




1 1 


n A ij ^n^- 


where 


Si = ^ 


Re[ Bij (k n ) 

I »[0 tj (k n ) ] 


and 


B = 
n 


1 0 -k 


Ok 0 


n 


2 1 2 
b + k 
1;j n 


b . .k 

1] n 

2 2 
b + k 

n 


2 2 

b . + k 

n JLj 3 n 


b .k 
n lj3 n 

2 , 2 
b . + k 

n ij 3 n 


(A6) 


Define the square of the magnitude of the error over all data points included in the 
least-squares fit to be 


nkk 


11 


Z T 

e e 
n . . n . 


(A7) 


n=1 


il il 


where nkk. < n,. Then the coefficients which provide the least-squares error sub- 

. IK ~ 

lect to the imposed constraints are found by determining the A^^ and >qj, which 
minimizes 


^ C ij 


= e. 


"il 
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(A8) 
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Performing the indicated minimization yields 


nkk- i 

v b t b ; 
/ > n n| 

— 

C ij 

n=1 • 


:"4r\ 

o 
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nkk. 
n=1 


1 3 


°ij 


k. 


(A9) 


Selectable Constraints 

The user may select constraints to impose upon the p-plane coefficients for 
column j from the following list: 

1. Force agreement with tabular value at k = 0 (asumes that k 1 =0): 


Oij (0) = Q ± j (0) = A 0i = [1 


0 0 




(all 1) 


2. Constrain slope at k = 0 based upon tabular data for column j (assumes 


that k^ = 0 and k 2 is small): 


0Q. 


13 


5p 


(P) 


p=0 


0 1 0 




Jr v 


Im [ P (k )] 

. = p 1— (all i) 

x 3 k o 


3. Constrain slope at k = 0 for column j to be related to the tabular data 

- R e[0 (k,- 0 )] 


in the mth column at k = 0 (assumes that k^ = 0) : 


dp. . 

13 

dp 


(p) 


p=0 


0 1 0 


’1j 


’■V 


A i j - 


lm 


This constraint can be used to relate plunge and pitch columns at k = 0 
(ref. 23) 


4. Constrain slope at k = 0 to be zero: 


dp 


-(P) 


p=0 


0 1 0 


( 1j 


n *j ] 


|Aij = 0 
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5. Null out A, : 

K j 

For example, let k = 2, 

A 0 = [0 0 1 ... 0] A . . = 0 

2 ij « 

where 1 occurs in column k + 1. 


6. Require a precise fit at a specific, nonzero frequency k which need not 
be a tabular point: 


Re[Q ij (k)] = 


1 0 -k 


,2 , 2 

b . + k 


2 , 2 
b n i + k 


A li = Re [Q. . (k) 


Im[Q i .(k)] = 


.0 k 0 


b lj k 


2 , 2 
b^ . + k 

id 


bn ljD k 


b n i + k 


A-- = Im[Q. . (k) ] 
1 J 1 J 


If this constraint is selected, the tabular data are interpolated to 
obtain Q^j(k). 


Instructions 
are given in 


for preparing program input to include one or more of these 
table II. See specifically the variables ICOF, THETAN, and 


constraints 

SPKO. 
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EXAMPLE OF INCLUDING SCHEDULED COMPONENTS IN COMPENSATOR 


STABCAR contains a provision for including compensator elements that vary as 
functions of Mach number and dynamic pressure. This capability is provided in func- 
tion subprogram SCHEDUL where control/sensor pair identification, Mach number, and 
dynamic pressure are available. A phase error <|k ^ introduced into the compensation 
for one control/sensor pair is also available. There are block data statements in 
SCHEDUL. The user inputs data required to define his scheduling laws by using the 
block data statements. He then provides code in SCHEDUL which utilizes the data and 
the previously mentioned parameters to define the desired scheduling laws. An exam- 
ple illustrating the procedure for a single-input/single-output control law follows. 
The scheduling law implemented is that shown in figure 7. 
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COMPLEX FUNCTION SCHEDUL < I C , JC , XMACH , QEAR , S ) 

COMMON /GAIN/ I GA I N , DE L G A I N ( 10,20) , GN( 10 , 20) , IPH,JFH, PHASE , PH ERR 
, ISCHDUL (10,20) 

COMPLEX PHERR 

COMMON/ AEROP /NMSAVE , NCSAVE , NRSAVE 
C , NKSAVE , NMNCS 

COMMON / GNSCHC / C OF 1 < 1 0 , 2 0 ) , CO F 2 ( 1 0 , 2 0 ) , C OF 3 ( 1 0 , 2 0 ) , COF 4 ( 1 0 , 2 0 ) 
,COF5< 10, 20) , COF 6 ( 10,20) ,COF7<10,20) ,COF8(10,20) ,COF9( 10,20) 

, COF 10(10,20) 

COMPLEX S 

COMMON/ SELECT/NMODES , NOC ( 25 ) , NKNEW , NCNEW , NRNEW 
C , NOK ( 15) , ISELECT 


DATA 

COF 1/0 . f 

0 . , 

61 , 197*0/ 

C 

, C O F 2 / 0 . 

, 0 . 

2 . 3924807E+04 , 1 9 7 *0 . / 

C 

, C O F 3 / 0 . 

, 0 . 

3 . 1473 1 9 E - G 5 , 197*0 . / 

C 

, COF4 / 0 . 

, 0 

,1034.21,197*0./ 

C 

, C O F 5 / 0 . 

, o . 

,1,197*0/ 

C 

, C O F 6 / 0 . 

, 0 . 

,735,197*0./ 

C 

, C OF 7 / 0 

, 0 . 

,8 3 3 . 3, 1 9 7 * 0 . / 

C 

, CO F 8 / 0 . 

, 0 . , 

, . 00 3 3 2 5 7 1 5 4,1 9 7 * 0 / 

C 

, C O F 9 / 0 . 

, 0 . . 

, 300 ,1 97 * 0 / 

C 

, COF 10/0 

, 0 

, , 75 . , 1 9 7 * 0 . / 


NMNC = NMODE S -NCNEW 
C F 1 = CO F 1 ( I C , JC ) 

CF2 = COF2 ( IC , JC > 

CF3=COF3 ( I C , JC ) 

CF 4 = COF 4 ( I C , JC ) 

CF5=COF5 ( IC , JC) 

CF6=COF6 ( IC, JC) 

C F 7 = COF 7 ( IC, JC) 

CF8 = COF8 ( IC , JC ) 

CF9=COF9 ( IC , JC) 

CF 1 0 = COF 1 0 ( IC , JC ) 

C*--> PHASE ERROR CAN BE INTRODUCED BETWEEN CONTROL I PH AND SENSOR JPH 

C* E XP ( I * PHERR ) 

SCHEDUL = ( 1,0.) 

IF ( IC . EQ . I PH. AND . JC . EQ . JPH) SCHEDUL =SCHEDUL * PHERR 
IF( I S CHDUL ( IC, JC) . NE . 1 ) RETURN 

C*--> THIS CODE IMPLEMENTS THE SCHEDULED COMPONENTS OF THE CONTROL LAW 

C* DEFINED IN FIGURE 7 

XM= XMACH 

C*--> I F ( M ( 0.61) THEN M= . 6 1 

1 F ( XM . LT . CF 1 > XM= CF 1 
QB=QBAR 

C*--> I F ( Q < 2.3924+04) THEN Q = 2.3924+04 

I F ( QB LT CF2 ) QB=CF2 

C*--> -1 (= KS = (3 . 1417-05 )*( 1 . 61-M) *(Q - 1 0 34 .2 1 ) (= 1 

X KS = C F 3 * ( 1+CF1-XM) * (QB-CF4 ) 

IF (XKS GT. CF5)XKS=CF5 

C*--> -75 (= DS = 735 . -833 . 33*M + 0.0033267*Q (= 300. 

DS = CF 6-CF7*XM+CF8*QB 
I F ( DS GT . CF9 ) DS=CF9 
I F ( DS LT . CF 1 0 ) DS=CF 1 0 

C*--> (KS/ (S+DS) ) 

SCHEDUL=CMFLX ( XKS , 0 . ) / ( S+DS ) ‘SCHEDUL 

RETURN 

END 
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OPTIONS FOR HANDLING DENOMINATOR DYNAMICS ARISING 
FROM ACTUATORS , FILTERS , AND SENSORS 

This appendix describes options contained in STABCAR which allow denominator 
dynamics arising from actuators, filters, and sensors to be cleared from the denomi- 
nators of the equations. This option is employed when control roots are to be found 
with gains near zero to avoid division by zero. 

Consider the following transfer function between actuator i and sensor j (see 
eqs . (9) and (13) ) : 


(6 >. 



r . 

D 


\f-i <t>. 


= (T ) . . G. . (T ) . . e 
A ll lj LS lj 


1 1 


( T ' ) . . ( T ) , . 
L i 3 S 


Let 


3 i j - G i j (t LS h j e 






and 


T ij ~ ( t A ^ ii ^ t L ^ i j ^ T S ) j j 


(Cl ) 


(C2 ) 


Define 


d-d, 


d^ s 


d 2 s2 + 


+ d 


(C3 ) 


as the common denominator of 



n$, j — 1, ..., Dg . De f i ne 


1 D 


= LjLj d 


(C4 ) 


An option is provided to multiply the control rows of either equation (1) or (2) 
through by d. The multiplication is performed if the variable IDMULT is input as 1. 
If IDMULT is input as 0, the multiplication is not performed. If equation (1) is to 
be employed, the variable ICSACT is input as 0. If equation (2) is to be employed, 
ICS ACT is input as 1. 

When characteristic roots are determined with the determinant iteration 
approach, the modifications to equations (1) and (2) that occur for IDMULT - 1 are 


28 



APPENDIX C 


evident. The modifications employed for the matrix iteration option are defined 
below (d r is the first nonzero coefficient in the polynomial d(s) (see eq. (C3) ) 

1. Let ICSACT = 1. Partition the matrix (see eqs. (2) and (18)) 


A = 




A 6H, A 66 


(C5 ) 


A^^ and A^ have the following definitions: 


(a) For IDMULT = 0, 


. 2 N ( s ) H 

*65 " 5 dTFT 


*65 - 0 


(b) For IDMULT = 1, 




s N(s)H 

d s r 
r 


-s 2 [d(s) - d s r ] 


66 


d s 
r 


2. Let ICSACT = 0 (see eqs. (1) and (18)). 
(a) For IDMULT = 0, 


A = - 


M. . 
11 


/V M - pM..l\s 2 + lg.(s) K. 

[I L 1 


+ K + qQ (n ,— ) 

\ rla U j 


- F, 


D 


N (s ) H 

dur + t he (s) h e 


(b) For IDMULT = 1, 


A = 


-1 


d s 
r 


li 


M - 


- f d t he < s) V 


M. . 
11 


s + g. 


c 

g. (s ) k. . 
i n 


+ k + qQ 


/ sb 
? Ka'U 

(d (s ) - d s") 


d(s) - F d N(s) HV - s 


d s 
r 
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ALGORITHM FOR MATRIX ITERATION FOR CHARACTERISTIC ROOTS 

The matrix iteration algorithm for characteristic root determination employed in 
STABCAR during convergence difficulties and for obtaining initial estimates for the 
characteristic roots is presented in this appendix. 


Let 


0 

P 


n \ 


option type for which matrix iteration routine is being employed. 
0 = 1 during convergence difficulties; 0^ = 2 if initial 

estimates are being sought and all ^ - 0; Op = 3 if initial 

estimates are being sought and some 


number of characteristic roots to be found by matrix iteration. n^ = 1 
if Op = 1 ; n^ = n^ if 0 p = 2 and ICSACT = 0; n^ = n^ + n fi 

if 0 = 2 and ICSACT = 1 ; n. = Number of nonzero s-; if 0 = 3 
p A. x 2 P 

counter for number of roots being found. i = 1, . .., n^ 


i root number in determinant iteration for which convergence difficulties 

n 

were incurred 


j counter determining number of iterations for a root 

JL eigenvalue number (from set of eigenvalues of which 

corresponds to , the initial estimate of current root being 

determined (see step 2) 


' 

11 ^ difficulties were incurred 

ith root being determined 

j th estimate for s^ 

square root of eigenvalue of &(s' 
i.e. , s' ± = 

j 


s - 
1 


GO 


best initial estimate for ith characteristic root 
second initial estimate for ith root (see step 6) 

estimate from determinant iteration routine for root for which convergence 


which corresponds to s^ , 


frequency used in estimating s^j if Op = 2 
see step 2 
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Step 1 


Step 


Step 3 


Step 


Step 


inverse of generalized mass matrix 
square root with positive imaginary part 
Initialize: 


i = 1 


s- = S/. x if 0 ~ 1 

X 1 (l n } 2 P 


= (-0.5to , 03 ) if 0=2 
r r p 


= s- if 0=3 
1 2 p 

Determine 


Se t j = 1 : 

Determine eigenvalues {\^} of ^ (see eqs. (18)) 

Determine index 1 such that s ^ 


£ = n y - (i - 1 ) if 0 = 2 

a P 


Or 1 is defined such that s! is closest in magnitude to s-' if 

1 1 


°p * 2 


: Determine convergence: 




< 0,0001 


: Set 




if 


j 


1 


go to step 6 (convergence) 




* a i 



-1 


a » i 



if (j > 1) 


: Set j = j + 1 


If j > 30, print out nonconvergence, and go to step 7 
Determine eigenvalues {A^} of (see eqs, (18)) 

Set s! = 

j 

Repeat step 3 
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Step 6: 

Set 

S <V2 = "ij 

and return if 0 = 

P 


Set 

Si 2 = s[ . and 

Si 1 = 0.99s-j_ 2 if 

Step 7: 

Set 

i = i + 1 



If 

i > return 



Set 

s’, = \| X U-1) 

if O p = 2 


* 

s il 

= s \ 2 ° P = 

3 


Repeat step 2 
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GRAPHICS ROUTINES NOT INCLUDED IN STABCAR 

This appendix is devoted to the description of elements of the Langley Graphics 
package that have been employed by STABCAR but are not included as a part of STABCAR. 
The description is included to facilitate a substitution of equivalent routines. 
Alternatively, the Langley Graphics package can be purchased separately from Computer 
Software Management Information Center (COSMIC), Suite 112, Barrow Hall, University 
of Georgia, Athens, GA 30602, as LAR 12722 entitled LRCGOS (Langley Research Center 
Graphics Output System). 


SUBROUTINE PSEUDO 

PURPOSE: To create and write an appropriately named Plot Vector File. Through link- 

ages set up by an initial call to PSEUDO, all subsequent graphics data generated by 
the user will be routed through one of the PSEUDO entry points and written on the 
Plot Vector File. The PSEUDO processor is designed for use with the frame depen- 
dent postprocessors described in Section 1.3, Volume IV, of the Computer Programing 
Manual. 

USE : CALL PSEUDO 

This will establish a Plot Vector File named SAVPLT. 

CALL PSEUDO (6LMYFILE) 

This will establish a Plot Vector File named MYFILE. 

NOTE: The Plot Vector File (or Files) will usually be written to disk (as 

opposed to tape) and may be postprocessed following user program termination 
via appropriate specification of one or more PLOT control cards. (See 
Section 1.3, Volume IV, Computer Programing Manual). 

RESTRICTIONS : An initializing call to PSEUDO (with or without a file name argument) 

must be made prior to any calls to CALPLT or any other graphics output routine. 

OTHER CODING INFORMATION: CALPLT and NFRAME are entry points in PSEUDO. PSEUDO 

checks for indefinite and infinite value and provides a traceback. 


ENTRY NFRAME 

PURPOSE : Tb provide users specific means of executing frame advance movements on any 

plotter device via an appropriate frame-oriented device postprocessor. Frame 
advance distances are generally defined to be incremental from current frame ori- 
gin. A frame is defined as a finitely bounded plotting area? it is an imaginary 
boundary around the usable plotting area. CALL NFRAME is intended to be used as a 
frame advance mechanism, not as a plot origin offset. 

USE: CALL NFRAME 
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ENTRY CALPLT 

PURPOSE: To move the plotter pen to a new location with pen up or down. This is the 

basic plotting routine that is used by all other graphics routines; therefore,- a 
check is made on the parameters for indefinite and infinite. If any of the parame- 
ters are indefinite or infinite, the program will abort and a traceback to the 
routine is provided. (See "OTHER CODING INFORMATION.") 

USE : CALL CALPLT ( X, Y, I PEN ) 

where 

X, Y are the floating-point values for pen movement. 

IPEN = 2 pen down 

= 3 pen up 

Negative IPEN will assign X = 0, Y = 0 as the location of the pen 
after moving the X,Y (create a new reference point). 

RESTRICTIONS: All X and Y coordinates must be expressed as floating-point inches 

(actual page dimensions) in deflection from the origin. 


SUBROUTINE LEROY/BALLPT 

PURPOSE: To set up the parameters necessary to accommodate plotting with the liquid 

ink pen. Once set, this mode will remain. 

USE: CALL LEROY 

RESTRICTIONS : The CALL LEROY is only recognized by the CalComp postprocessor; it is 

ignored by the other postprocessors. In addition to reducing the speed of the 
plotter for all plotting movements, the number of plot vectors in any annotation is 
considerably increased. 

The CALL LEROY must be made prior to any plotting calls, but after 
the CALL PSEUDO. 
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SUBROUTINE NUMBER 

PURPOSE: To convert a floating-point number to BCD (expressed in F format) and to 

draw the resulting alphanumeric characters. 

USE: CALL NUMBER ( X, Y, HEIGHT, FPN, THETA , NODIGIT) 

where 

X,Y are the coordinates in floating-point inches of the left lower corner of the 
first digit of output. 

HEIGHT is the height of the plotted number in floating-point inches. (See NOTATE 
routine. ) 

FPN is the floating-point number to be drawn. 

THETA is the angle in floating-point degrees at which the number is to be drawn. 

(See NOTATE routine. ) 

NODIGIT is the number of decimal digits to the right of the decimal point for 
output. 

N0DIGIT=-1 or NODIGIT=0 both specify no decimal places; however, -1 sup- 
presses the decimal point. 
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SUBROUTINE NOTATE 

PURPOSE: To draw alphanumeric information for annotation and labeling and provide 

special centered symbols for annotation of data points. 

USE: CALL NOTATE ( X, Y, HEIGHT, BCD, THETA, NOCHAR) 

where 

X, Y are the floating-point coordinates of the first character. 

For alphanumeric characters, the coordinates of the lower left-hand 
corner of the characters are specified. 

HEIGHT specifies character size and spacing in floating-point inches for a 

full-size character. 

BCD is the string of characters to be drawn and is usually written in the 

form: nHXXXX (the same way an alpha message is written using 

FORTRAN format statements). Instead of specifying alpha information 
as above, the beginning storage location of an array containing alpha- 
numeric information may be given. 

THETA is the angle in floating-point degrees at which the information is to be 

drawn. Zero degrees will print horizontally reading from left to 
right, 90° will print the line vertically reading from bottom to top, 
180° will print the line horizontally reading from right to left 
(i.e., upside down), and 270° will print vertically reading from top 
to bottom. 

NOCHAR is the number of characters, including blanks, in the label. NOCHAR 

is limited to 400 characters in one call. 
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SUBROUTINE AXES 

PURPOSE: To draw a line, to annotate the value of the variable at specified inter- 

vals with or without tic marks, and to provide an axis identification label. 

USE: CALL AXES ( X, Y, THETA, S ,ORG, SPX, TMAJ, TMIN, BCD, HEIGHT, NOCHAR ,NDEC ) 

where 

X, Y are the coordinates in floating-point inches of the starting point of 

the axis with reference to the plotting area origin as established by 
CALPLT. 

THETA is the angle of rotation measured counterclockwise from the X-axis in 

floating-point degrees. Normally, THETA is 0° for an X-axis and 90° 
for a Y-axis. 

S is the length of the axis in f loating-point inches. Should be a multi- 

ple of TMAJ. 

+S will generate tic marks. 

-S will eliminate tic marks. 

ORG is the functional value to be assigned to the origin (i.e., the value of 

the first scale) in floating point. 

SFX is the adjusted scale factor for the array to be plotted (change in 

value per inch). 

NOTE: Values of ORG and SFX which will produce a reasonable scale may 

be calculated using subroutine ASCALE or BSCALE. 

TMAJ is the distance in floating-point inches for major tic marks (0.25 inch 

high). If the values are integer multiples, the decimal point and 

decimal places are eliminated. A negative TMAJ will cause the actual 

value to be written instead of the adjusted value. 

TMIN is the number of divisions per inch in floating point for minor tic 

marks (0.125 inch high). To eliminate minor tic marks the following 
may be used: 

TMIN = 0. 

BCD is the character label for the axis (see NOTATE routine). 

HEIGHT is the height of the full-size characters in the BCD title. Numbers at 

the tic marks will be (0.75 ^HEIGHT) high. HEIGHT is in floating- 
point inches. 

If HEIGHT = 0. , all annotation will be eliminated. 
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NOCHAR is an integer specifying the number of characters in BCD title. A nega- 

tive NOCHAR places the annotation on the clockwise side of the axis 
and a positive NOCHAR places the annotation on the counterclockwise 
side of the axis, NOCHAR = 0 is not allowed. Lf it is desired to 
have no label, then the BCD parameter should be 1H, and NOCHAR either 
+ 1 or -1. Normally NOCHAR = 4 - for Y-axis and NOCHAR = - for X-axis. 

NDEC is the number of decimal places after the decimal point for the numbers 

placed on the axis at the major tic marks. NDEC may be omitted. If 
NDEC is omitted, the default value of 2 is set. 


ENTRY ASCALE 

PURPOSE: To compute a scaling factor for an array of numbers to be plotted over a 

certain area and find the minimum data value with the array. 

USE : CALL ASCALE( ARRAY, S ,N , K, DV) 

where 

ARRAY is the name of the array containing the floating-point values to be 

scaled. (See RESTRICTIONS.) 

S is the length (floating-point inches) over which the data are to be 

plotted (usually the length of one of the axes). 

N is the number of data values in ARRAY from which points are to be plot- 

ted in accordance with K. 

K is the interleave factor which specifies the sequence in which data are 

stored. 

= 1 indicates the values are stored sequentially. 

= 2 indicates that values are stored in every other location in the 
array, etc. 

DV is the number of divisions per inch of the plotting paper to be used 

(should be: 10.0 or 20.0). 

RESTRICTIONS : The array must be dimensioned to include storage space for two extra 

elements per interleave factor. For example: N = 100, K = 1, DIMENSION 

ARRAY (102); N = 75, K = 3, DIMENSION ARRAY (231). 

METHOD: This routine scans the elements starting at location L of the array to find 
the minimum and maximum. (If the scaling begins with first location of array, 

L - 1 ; if scaling begins with third location of array, L = 3. ) ASCALE computes an 
adjusted minimum (origin) and stores it in ARRAY ( L+ ( N*K) ) and computes a scale 
factor and stores it in ARRAY( L+( N*K) +K) • Hie scale factor will be (A*10**J) 
where A is 1, 2, 4, or 5 and J is integer power. The data in the array may be 
scaled to f loating-point inches by: 

SV = (AE - AMV ) /SF 
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where 

SV = scaled value 

AE = present value of array element 
A MV = adjusted minimum value computed by ASCALE 
SF = scale factor computed by ASCALE 


APPENDIX E 


SUBROUTINE BSCALE 

PURPOSE: To compute a scaling factor for an array of numbers to be plotted over a 

certain area and to find an acceptable minimum data value for the array. 

USE: 1. Stores the adjusted minimum value (origin) in ARRAY ( L+( N*K) ) and scale fac- 
tor in ARRAY ( L+ (N*K) +K) • If the scaling begins with first location of array, 

L = 1 ; if scaling begins with second location of array, L = 2, etc. 

CALL BS CAL E( ARRAY, S ,N, K, TMAJ, M,ORG) 

where 

ARRAY is the name of the array containing the floating-point values to be 

scaled. (See RESTRICTIONS.) 

S is the length (floating-point inches) over which the data are to be 

plotted (usually the length of one of the axes). 

N is the number of data values in ARRAY from which points are to be plot- 

ted in accordance with K. 

K is the interleave factor which specifies the sequence in which data are 

stored. 

= 1 indicates that values are stored sequentially. 

= 2 indicates that values are stored in every other location in the 
array, etc. 

TMAJ is the distance in floating-point inches for major tic marks along the 

axis (length of S). 

= 0. the routine scales in a manner similar to routine ASCALE. 

> 0. the routine will adjust the origin and scale factor in order to 
give regular multiples of the delta number that will be written at the 
tic marks by the AXES routine. 

M is an integer that has multiple uses. 

= 0 if M is zero and TMAJ is zero, the routine scales in a manner simi- 
lar to routine ASCALE. 

> 0 if M is positive and TMAJ is zero, the data will be scaled over the 
full range of S, with no attempt to adjust the scale factor for 
appearance. 

< 0 if M is negative, the valid value contained in parameter ORG will be 
used for a forced origin. 
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ORG is a floating-point value which will be used for a forced origin if 

parameter M is negative. The value in ORG may not be exactly the same 
on return from subroutine. 

ORG may be omitted if M is positive or zero. 

RESTRICTIONS : The array must be dimensioned to include storage space for two extra 

elements per interleave factor. For example: N = 100, K = 1, DIMENSION 

ARRAY (102); N = 75, K = 3, DIMENSION ARRAY (231). 
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SUBROUTINE PNTPLT 

PURPOSE: To draw NASA Standard Plot symbols centered on a given coordinate value. 

USE: CALT. PNTPLT(X, Y, ISYM, IS) 

where 

X is the X coordinate for the centered symbol in floating-point inches. 

Y is the Y coordinate for the centered symbol in floating-point inches. 

ISYM is an integer specifying the symbol to be used. (See figs. El and E2. ) 

= 21 for a point. 

= 22 for a plus sign +. 

IS is an integer value specifying the size symbol to be used. 

= 1 small 
= 2 me di um 
- 3 large 
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Figure El.- Plot symbols. 
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SUBROUTINE LINPLT 

PURPOSE; To draw a line between and/or draw a NASA standard symbol at each succes- 
sive data point (stored in an array), or to draw a continuous line through the data 
points and draw a specified symbol at the end of the line. 

USE : CALL LINPLT( X ARRAY, Y ARRAY, N, K, J, ISYM, IS,NG) 

where 

XARRAY, are the names of arrays containing the X values and Y values, respec- 
YARRAY tively, to be plotted. Values must be in floating point. (See 

RESTRICTIONS. ) 

N is the number of points to be plotted. 

K is the interleave factor which specifies the sequence in which data are 

stored. 

= 1 indicates that values are stored sequentially. 

= 2 indicates that values are stored in every other location in the 

array, etc. 

J is positive for line and symbol plot, negative for symbol-only plot. 

Ihe magnitude specifies the alternate number of data points at which 
to plot a symbol. 

= 0 for line plot with a specified symbol at the end of the line. 

= 1 for symbol for every data point. 

= 2 for symbol for every other data point, etc. 

IS YM is an integer describing symbol to be used, see PNTPLT routine for list 

(fig. El). 

= 0 no symbol will be drawn. 

IS is an integer specifying the size of the symbol (see PNTPLT routine). 

NG is an integer used to specify data points which will be enclosed by a 

1-inch square grid, 20 x 20 lines to the inch. NG may be omitted. If 
NG is omitted, the default value of 0 is set. 

= 0 no grids. 

= 1 for every data point. 

= 2 for every other data point, etc. 

RESTRICTIONS: LINPLT expects the adjusted minimums and scale factors as described in 

ASCALE. 
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SUBROUTINE DLINPLT 

PURPOSE: To draw a dashed line between and/or draw a NASA standard symbol at each 
successive data point (stored in an array), or to draw a continuous dashed line 
through the data points and draw a specified symbol at the end of the line. 

USE : CALL DLINPLT (XARRAY, Y ARRAY, N, K, J, ISYM, IS,LP) 

where 

XARRAY, are the names of arrays containing the X values and Y values, respec- 
YARRAY tively, to be plotted. Values must be in floating point. (See 

RESTRICTIONS. ) 

N is the number of points to be plotted. 

K is the interleave factor which specifies the sequence in which data are 

stored. 

= 1 indicates that values are stored sequentially. 

= 2 indicates that values are stored in every other location in the array, 

etc. 

J is positive for line and symbol plot, negative for symbol only plot. The 

magnitude specifies the alternate number of data points at which to plot a 
symbol. 

= 0 for line plot with a specified symbol at the end of the line. 

= 1 for symbol for every data point. 

= 2 for symbol for every other data point, etc. 

ISYM is an integer describing symbol to be used, see PNTPLT routine for list 
(fig. El). 

= 0 no symbol will be drawn. 

IS is an integer specifying the size of the symbol (see PNTPLT routine). 

LP is an integer used to specify the line pattern. 

= i 

= 2 

= 3 

= 6 
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= 7 
= 8 


RESTRICTIONS: DLINPLT expects the adjusted minimums and scale factors as described 


in ASCALE 
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DETAILED FUNCTIONAL DESCRIPTION OF OVERLAYS , MAJOR SUBROUTINES, 

AND IDENTIFICATION OF KEY PARAMETERS 

In this appendix, a description of the functions of each overlay and its major 
subroutines is given including definitions of key input and internally generated 
parameters that control program flow. 

Overlay MAIN 

Overlay MAIN is the STABCAR executive. Data are input here which determine the 
type of case to run, where to find array data, and which overlays are required. Key 
parameters are defined as follows but are defined more completely in table II: 

ICASE Integer indicating type of case (interactive or batch) to run and where 

to get data. 

ISPLANE Integer indicating whether to make p-plane approximation to aerodynamic 

forces. 

I SENSE Integer indicating whether to perform interpolation for modal 

deflection or slopes at specified sensor locations. 

INTERP Integer indicating whether to compute geometric coefficients required 

for surface spline interpolation. 

IPLATE(I) Integer defining whether Ith structural section is modeled as beam or 

as plate. 

ICSREAD Integer defining whether sensor data are to be read from TAPE5. 

ICONSYS Integer defining whether control system is to be included. 

ICHSEN, ICHFIL, Integers defining whether changes to previously defined control 

and ICHACT system are to be made in current case; if not, overlay CONTROL 

is not needed. 

I SELECT Integer defining whether subset of current model is to be selected and 

where data for current model can be found. 

IFLUT Integer defining whether to perform stability analysis in current case. 

IAPLT Integer defining whether to make plots of aerodynamic forces. 

IPKPLT Integer defining whether to save data for making plots of stability 

characteristics. 

IDELPLT Integer defining whether to delete unwanted characteristic root data 

from data storage file TAPE1 . 
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Overlay INTERP 

Geometry-related, mode -independent coefficients which are needed for a surface 
spline fit for later interpolation purposes in overlay SENSOR for plate-type sections 
are computed by overlay INTERP. Primary subroutines called are INTEXT, SYMSURF, 
NONSURF, and MATINV. Key parameters are defined as follows: 


NNODES (or NNODS) Total number of nodes at which modal data are available. 
NPMX Maximum number of nodes in one structural section. 

NSECTNS Number of structural sections vehicle is divided into. 


I PLATE 
IPLATE (I) 


NSECTNS x 1 array. 

1 Ith section is modeled as a plate. 
0 Ith section is modeled as a beam. 


TAB NNODS x 4 array which, for plate-type sections, defines x, y coordinates 

(global) and z-def lection; fourth column is not employed for plate 
sections. For beam sections, array elements are distances along axis 
of rotation at which data are available, z-def lection, rotation about 
swept Y-axis (y* ) and rotation about swept X-axis (x* )• (See 
fig. 3. ) There is a TAB array for each elastic mode shape. 

ISS 2 x NSECTNS array. ISS(1,I) is node number of first node of Ith 

structural section; ISS(2,I) is node number of last node of Ith 
structural section. 


RO, XO, and YO NSECTN x 1 arrays. For plate-type modes, RO = 1/B, where B is 

one-half root chord of Ith structural section; XO(I) is 
x-coordinate of Ith section root semichord; and YO(I) is 
y ^coordinate at root semi chord. For beam-type sections, 

IPLATE = 0, (X0( I) , Y0( I) ) is elastic axis intercept with Ith 

structural section root chord; RO is sweep angle of elastic 
axis of Ith structural section in radians. 

Array of geometry coef f icients. 

Subroutine INTEXT .- Subroutine INTEXT is an extension of INTERP and it is here 
that the TAB arrays are actually read and geometry coefficients are obtained. Sub- 
routines called are SYMSURF and NONSURF. Key parameters used are all those pre- 
viously defined for overlay INTERP. 

Subroutines SYMSURF and NONSURF .- Subroutines SYMSURF and NONSURF compute the 
mode-independent geometry coefficients for symmetric and antisymmetric or nonsym- 
metric plate sections, respectively. Subroutine called is MATINV. Key parameters 
used are TAB and H. 


t 


Internally 


generated parameters. 
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Subroutine MATINV .- Subroutine MATINV performs the matrix inversion needed in 
computing H, the geometry coefficients. 


Overlay SENSOR 

Overlay SENSOR determines surface spline coefficients for fitting mode shapes 
over plate-type sections. Beam-type modal data are also input here. Interpolation 
of beam and plate modal data is performed to obtain modal coefficients defining 
either linear or angular displacements that would be measured by sensors at specified 
locations. These modal coefficients are stored for subsequent use and define the 
H matrix in equation (1) or (2). Primary subroutines are SENCTL, DEFLECT, CSIUNI, 
IUNI, SPLDER, ZFUN2, AFUN2, ZFUN, AFUN, ZDER, and ADER. Key parameters are defined 
as follows: 


NS 

XS and YS 
ITYPE 
NSS 
ISS 

IPLATE ( I ) 
TAB, R0,X0, 


Number of sensors. 

NS x 1 array of sensor locations. 

f = 1 Linear measurement. 

- 2 Angular measurement. 

NS x 1 array identifying structural section where each sensor is 
located. 

2 x NSECTNS array defining beginning and ending node numbers for each 
structural section; see section "Overlay INTERP" for additional 
information. 

1 Ith structural section is modeled as plate. 

= 0 Ith structural section is modeled as beam, 
and YO Defined in section "Overlay INTERP." 


NR 


XCG 


U and V 
CS 

w 


Number of rigid-body modes. In this overlay NR must be either 0, 1, 

or 2. If NR = 1 , the rigid-body mode must be the negative of 
plunge; the user may input his own rigid-body modes by setting 
NR = 0 and treating his rigid modes as if they were elastic. 

x-coordinate of airplane center of mass. Since first rigid-body mode 
is assumed to be the negative of plunge, and second, pitch about 
center of mass, this is only required if NR = 2. 

Sensor location with respect to swept axis. 

Computed n g x n^. array containing modal participation of each 
particular mode to input of each sensor. 

Geometry coefficients computed in overlay INTERP for plate-type 
sections. 
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Subroutine SENCTL .- In subroutine SENCTL, the data required for modal interpo- 
lation are read in from file TAPE1 0, and after linear or angular displacements are 
calculated, resulting modal coefficients are written to data storage file TAPE1 . 
Primary subroutine called is DEFLECT. Key parameters used are all those defined for 
overlay SENSOR. 

Subroutine DEFLECT .- In subroutine DEFLECT, the actual modal coefficients defin- 
ing linear or angular displacement inputs to sensor are obtained. Primary subrou- 
tines called are CSIUNI, I UNI, SPLDER, ZFUN2, AFUN2, ZFUN, AFUN, ZDER, and ADER. Key 
parameters used are all those defined for SENSOR except NS and NSS. 

Subroutines IUNI and CSIUNI .- Subroutines IUNI and CSIUNI are used to perform 
interpolation of modal data used in calculating both linear and angular displacements 
for beam-type sections ( I PLATE = 0). IUNI performs linear and quadratic interpola- 
tion when there are less than four tabular values and CSIUNI performs a cubic spline 
interpolation when there are four or more tabular values. Key parameters used are 
TAB and V. 

Subroutine SPLDER . - . Subroutine SPLDER interpolates for Ry* and estimates 
dRy'/dy* for beam-type sections ( I PLATE = 0) used in computing angular displacements 
( ITYPE = 2). Key parameters used are TAB and V. 

Subroutines ZFUN 2 and AFUN 2 . - In subroutines ZFUN 2 and AFUN 2, spline coeffi- 
cients are determined for symmetric (ISYM = 1), nonsymmetric or antisymmetric 
( IS YM ^ 1) plate-type sections ( IPLATE = 1). Key parameters used are TAB and W. 

Subroutines ZFUN and AFUN.- In subroutines ZFUN and AFUN, modal coefficients 
for linear displacement-type sensors are determined for symmetric (ISYM = 1) and 
nonsymmetric or antisymmetric (ISYM ^ 1) plate-type sections (ITYPE = 1 and 
IPLATE = 1). Key parameters used are TAB, W, XS, and YS. 

Subroutines ZDER and ADER .- In subroutines ZDER and ADER, modal coefficients 
for angular displacement or rotational-type sensors are determined for symmetric 
(ISYM = 1) and nonsymmetric or antisymmetric (ISYM 1) plate-type sections 
(ITYPE = 2 and IPLATE = 1). Key parameters used are TAB, W, XS, YS, and R0. 


Overlay MATINPT 

Large data arrays such as the aerodynamic data, generalized masses, generalized 
stiffnesses, and control system dynamics definitions are input in overlay MATINPT. 

In addition, the selection feature is available which enables one to choose a mathe- 
matical model that is a subset of a previously defined model. Primary subroutines 
called are MAT EXD, CSINPUT, and SELECTS. Key parameters are defined as follows: 


ISPLANE 



1,2 

0,3 


Read p-plane coefficients. 

Read aerodynamic coefficients. 
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I SELECT 


ICONSYS 


KASE + 


ICSREAD 


INTERAC 


ICHSEN 


ICHFIL 


ICHACT 


NM and 
NMODES 


< 0 Read modal data from data storage file TAPE1 (select from 
previously used data). 

> 0 Read modal data as formatted input from file TAPE5 (select from 
original data). 

= 0 No selection of array data; for initial case, data from file TAPE5 
are used? in subsequent cases, array data from preceding case are 
used for current case. 

= 1 Call subroutine CSINPUT to input sensor and control system data. 


V^= 0 No control system input. 

Number of cases run during current program execution. 

f- 0 Read sensor data from data storage file TAP El , record 16 as 
J computed and stored in overlay sensor. 

1 Read sensor data as input from file TAPE5 • 

r - 1 Interactive run is being made (additional input after initial 
ization case is requested from keyboard). 






0 Run is in batch mode (additional input after initialization case 
must also be on TAPE2. 

1 Change in sensor dynamics. 

0 No change in sensor dynamics. 

0 No change in filter dynamics. 

1 Change in filter numerator. 


= 2 Change in filter denominator. 

3 Change in both numerator and denominator. 
r = 0 No change in actuator dynamics. 

= 1 Change in actuator numerator. 

= 2 Change in actuator denominator. 

= 3 Change in both numerator and denominator. 

Number of modes employed in model before selection and number to be 
selected ( ISELECT f 0). 


t 


Internally 


generated parameter. 


50 



APPENDIX F 


NR and Number of rigid-body modes employed in model before selection and 

NRNEW number to be selected ( ISELECT f 0)* 

NC and Number of controls employed in model before selection and number to be 

NCNEW selected (ISELECT f 0). 

NK and Number of reduced frequencies at which aerodynamic tabular values 

NKNEW are employed in model before selection and number to be selected 

(ISELECT ^ 0). If NKNEW = 0 (default), all will be used and no 
selection per frequency will be performed, 

NOC Array of mode numbers indicating which modes to select (ISELECT f 0). 

NOK Array of frequency indices indicating which reduced frequency data to 

select (ISELECT t 0 and NKNEW £ 0). 

NS Number of sensors . 

Subroutine MATEXD . - , A1 1 noncontrol array data are read in subroutine MATEXD and 
all array data are rewritten to data storage file, TAPE1 , in form to be used by other 
overlays. Primary subroutines called are CSINPUT and SELECTS. Key parameters used 
are ISPLANE, ISELECT, ICONSYS, KASE, INTERAC, NM, NR, NC , NS, and all arrays defining 
aeroelastic model and control system. 

Subroutine CSINPUT .- All control system data are read in subroutine CSINPUT, or 
default values set, including sensor deflection coefficients. Key parameters used 
are ICSREAD, ICHSEN, ICHFIL, ICHACT, NS, NC, and arrays defining control system. 

Subroutine SELECTS. - All modal and frequency selection of data is performed in 
subroutine SELECTS. Key parameters used are ISPLANE, INTERAC, KASE, NM, NMODES, NR, 
NRNEW, NC, NCNEW, NK, NKNEW, NOC, NOK, NS, and all arrays defining basic aeroelastic 
model. 


Overlay P PLANE 

The p-plane curve approximations to the oscillatory aerodynamic forces are com- 
puted in overlay PPLANE corresponding to a specified set of denominator constants 
C {b j } (eq. (3))). The best coefficients to employ in equation (3) for a given set of 
oscillatory aerodynamic force data at a number of reduced frequencies are determined 
in a least-squares sense. The resulting coefficients and the corresponding denomi- 
nator constants are stored for subsequent use in overlays PKFLUT and AEROPLT. Pri- 
mary subroutines called are SET and AERCOEF. Key parameters are defined as follows: 


NK 

NKK ( J) 


NCOEF 

NPOLYC(J) 


Number of reduced frequencies at which aerodynamic data are available. 

Jth element of NM x 1 array indicating that only first NKK reduced 
frequencies of aerodynamic data are to be employed in least-squares 
fit for Jth column of aerodynamic forces. 

Number of coefficients desired in p-plane fit, including polynomial and 
denominator (lag) terms (Max = 10). 

The number of polynomial coefficients for column J. 
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Integer indicating whether to include nth constraint (as defined below) when 
fitting Jth column of aerodynamic forces. (See appendix A.) It is assumed that 
= 0 and k 2 is small. 


ICOF(N, J) 



0 (Default) Nth constraint not employed in Jth column of 
aerodynamic forces. 


0 Employ Nth constraint. 


IC0F(1 , J) 
ICOF ( 2 , J ) 


> 0 Force agreement with first tabular value at = 0; 

Q. . (0) = Q. . (0) . 

ID ID 

> 0 For k^ = 0, force slope, 



dp 


p=0 


Im {Q l j ( e ) }/e , where e -^2* 


ICOF ( 3 , J ) 
, ICOF ( 4 , J) 

IC0F( 5, J) 

I 

IC0F( 6, J ) 
BN ( I , J ) 


> 0 For k^ = 0, Force slope = -Re ( ^ ) }/ ( b0^ ) . 

> 0 For = 0, Force slope = 0. 

= JL + 1 Force A^ = 0. 

j 

> 0 Force Qj(kg) = Q(kg) f° r some k Q , where Qj (kg) i- s interpolated 

value of Jth aerodynamic force at kg, obtained with a piecewise 

quadratic fit to Q.(k. ), i = 1, ..., n v . 

D 1 K 

Constant in Ith denominator term of form P/(BN + P) corresponding to 
Jth column (pressure mode)? if same set of BN is used for several 
columns they need only be input once. (See NCOL. ) 


THETAN Normalizing factor employed when applying constraint number 3, 

ICOF ( 3 , J ) = 1? Default = 1. 


SPK0 

! 


Value of reduced frequency for which p-plane approximation must 
precisely fit data. 


NM 


Number of modes being employed in model. 


NC 

NOC 


NCOL ( J) 


Number of controls being employed in model. 

Array of mode numbers indicating which modes are included in current 
mode 1. 

( - N Integer indicating preceding column number, if any, for which 
J BN ( I , J) = BN ( I , N) (default is 1). 

1^= 0 BN terms for Jth column are to be input. 
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NBN(J) f Number of nonzero BN for column J. 

C(I, J,K) Ith p-plane coefficient for (J,K) element of generalized aerodynamic 

force matrix. 

( = 1 Write percent errors for p-plane fit at tabular values to TAPE6. 

= 0 Do not write. 

Subroutine SET . - The number of lag terms is determined for each column of the 
aerodynamic forces based upon input parameters. The logic in subroutine SET was 
developed to reduce the input required of the user. The parameter NCOL(J) is 
employed to indicate whether lag terms for column J were read in (NCOL(J) = 0) or are 
the same as those for column JJ = NCOL(J). Key parameters used are NCOEF, BN, NM, 
NC, NOC, NCOL, and NBN. 

Subroutine AERCOEF .- Subroutine AERCOEF solves for the coefficients that yield 
the minimum error in the least-squares sense in the p-plane fit to the aerodynamics. 

A library subroutine MATOPS is employed to solve the set of linear equations which 
result. (See appendix A.) Key parameters used are all for overlay PPLANE. 


Overlay CONTROL 

For overlay CONTROL, a matrix of transfer functions T^T^T s is constructed. 

(See eq. (9).) The ( )— element relates the output of actuator i to the input 
of sensor j, excluding (G)^j and any scheduling or phase error. The polynomial 
coefficients for sensor, condensation, and actuator dynamics are constructed sepa- 
rately and then combined and a common denominator is found. Each element of the 
resulting matrix can then be multiplied by a distinct feedback gain and scheduling 
component in overlay PKFLUT. Primary subroutines called are CONTEXT, SENDYN, CLOGIC, 
FILDYN, ACTDYN, PMULT, NCOF, and EQUATEP. Key parameters are defined as follows: 


NC 

NS 

ISDYN(I) 


IORD(I) 


XKS ( J, I ) 


Number of controls. 

Number of sensors. 

^ = 0 Perfect sensor (default). 

L = 1 Sensor dynamics included. 
r = 0 Position sensor, 
v = 1 Rate sensor. 

2 Acceleration sensor. 

Jth coefficient of numerator polynomial for Ith sensor. 


t 


Internally generated parameter. 
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ASN ( J, I ) f 


ASD( J, I) 


r = 1 If J = IORD(I) + 1 and ISDYN(I) = 0. 

= 0 If J f IORD(I) + 1 and ISDYN(I) = 0, 

[^= XKS ( J J, I ) where J = JJ + IORD(I) + 1 if ISDYN(I) = 1. 
r 

= 1 If ISDYN(I) = 0 and J = 1 . 

/= 0 If ISDYN(I) = 0 and J > 1. 

= Jth coefficient of denominator polynomial for Ith sensor if 
L ISDYN(I) = 1. 


ISN(I) and Number of coefficients in numerator and denominator, respectively, 
ISD(I) of transfer function modeling Ith sensor dynamics. 


IACT(I) 


"= 0 Ith actuator is modeled as perfect actuator. 

,= 1 Actuator dynamics are included for Ith actuator. 


AACTN ( J, I ) and Jth coefficient of numerator and denominator of polynomial, 

AACTD( J, I ) respectively, of transfer function used to model Ith actuator. 


IAN ( I ) and Number of coefficients in numerator and denominator of Ith actuator 
I AD ( I ) transfer function. 


ACLN(K,I,J) and Kth coefficient of numerator and denominator polynomial, respec- 
ACLD ( K, I , J) tively, of transfer function which relates Ith actuator input 

to Jth sensor output excluding gain, scheduling, and phase 
error contributions. 


ILN ( I , J ) and Number of numerator and denominator coefficients in transfer func- 

ILD ( I , J) tion relating Ith actuator input to Jth sensor output. 


IFILTER^ j (K) Indicates which types of filters are to be combined to define 

overall filter transfer function relating actuator \ input to 
sensor j output. Maximum of 12 filters for each control- 
sensor pair. I FILTER ( K) is Kth filter type (sequence must end 
with 1 ) . 

S' 

- 1 No filtering. 

= 2 Notch filter. 

= 3 Integral filter. 

IFILTER ( K) / 

= 4 Proportional plus derivative filter. 

= 5 Lead-lag (or lag-lead) filter. 

= 6 General filter comprised of numerator and denominator 
v polynomials. 


t 


Internally generated parameter. 
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WNI(K) and Natural frequency and damping ratio of notch filter if 
SNI(K) IFILTER ( K) = 2. 

KIBDO(K) Gain employed in integral filter if IFILTER(K) = 3. 

KIBD1 (K) and Gains employed for proportional plus derivative filters, respec- 
KIBD2(K) tively, if IFILTER(K) = 4. 

AN 1 ( K ) and Time constants of zero and pole in lead-lag (lag-lead) compensator 

ADI ( K) if IFILTER(K) = 5. 

AFN ( I , K) and Numerator and denominator coefficients (constant term first) of 

AFD( I , K) general-type transfer function if IFILTER ( K) - 6 . 

ACTLN(K, I, J) + and Kth confuted coefficient of numerator and denominator polynomial, 
ACTLD ( K, I , J) ^ respectively, of overall transfer function ( T^TlT s ) p j 

relating actuator I output to sensor J input, excluding 
gain, scheduling, and phase error contributions* 

ICLN( I, J) ^ and Number of coefficients in numerator and denominator of overall 

ICLD(I,J) + transfer function (T T* T ). . . 

A L s 13 

DCT + Array of coefficients in common denominator d of all transfer 

functions (T T’T ). ., i.e. , common multiple of all ACTLD arrays. 

A L s lj 

NCT Number of coefficients in DCT • 

CNTLN(K, I , J) + Kth coefficient of polynomial corresponding to product of 

ACTLN ( I , J) and common denominator DCT = (T T'T ). .d. 

A L s 13 

Subroutine CONTEXT. - In subroutine CONTEXT the overall matrix product of the 
sensor, control logic, and actuator matrices of transfer functions is formed and a 
common denominator is found. Primary subroutines called are SENDYN, CLOGIC, and 
ACTDYN. Key parameters used are all except IFILTER, WNI, SMI, KIBOO, KIBOI, KIB02, 
AN1, ADI, AFN, and AFD. 

Subroutine SENDYN .- Subroutine SENDYN determines the numerator and denominator 
polynomial coefficients and order of the transfer functions relating sensor output to 
sensor input. Ihe numerator coefficients are adjusted for the appropriate power of 
s to define a position, velocity, or acceleration sensor. Subroutine called is 
NCOF • Key parameters used are NS, ISDYN, IORD, XKS, ASN, ASD, ISN, and ISD. 

Subroutine CLOGIC .- Subroutine CLOGIC sets up the input for calls to FILDYN in 
order to obtain the numerator and denominator polynomial coefficients which make up 
the transfer functions relating sensor output to actuator input. Primary subroutine 
called is FILDYN. Key parameters used are ACLN, ACLD, ILN, and ILD. 

Subroutine FILDYN . - Subroutine FILDYN contains the filter types that may be 
combined to construct control logic dynamics. Ihe mathematical expressions that 
represent the various filter types are defined in the section "Control-System Repre- 
sentation." Any combination of the filter types may be utilized to generate an 
overall filter. For example, if IFILTER ( 1 ) = 3, IFILTER ( 2) = 6 , IFILTER ( 3 ) = 6 , 

* Internally generated parameter. 
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and IFILTER( 4 ) = 1, the overall filter is the product of one integral and two general 
polynomial filters. When IFILTER ( K) = 1 is encountered, no further filter dynamics 
are computed for the particular control sensor pair under consideration. Therefore, 
the IFILTER array must end with 1. Subroutines called are NCOF and PMULT. Key 
parameters used are IFILTER, WN1, SNI, KIBDO, KIBD1 , KIBD2, ANI, ADI, AFN, and AFD. 

Subroutine ACTDYN • - Subroutine ACTDYN determines the order of the polynomials in 
the transfer function representing the actuator dynamics relating actuator output to 
actuator input. Subroutine called is NCOF. Key parameters used are IACT, AACTN, 
AACTD, IAN, and I AD. 

Subroutine NCOF .- NCOF is a function subroutine which returns the number of 
coefficients in a polynomial; that is, it determines the order +1 of a polynomial. 

Subroutine PMULT .- The utility subroutine PMULT multiplies two polynomials 
together. All polynomial arrays are assumed to have constant terms first. 

Subroutine EQUATEP .- The utility subroutine EQUATEP sets one polynomial equal to 
another. 


Overlay PKFLUT 

The matrix in equation (1) or (2) is formed by overlay PKFLUT. A user-specif ie d 
number of characteristic roots are determined as a function of variation in velocity, 
altitude, density, or feedback gains and the resulting output is stored for possible 
subsequent plotting (if IPKPLT f 0) . (See section "Characteristic Root Determina- 
tion." ) By operating in a multicase mode, the effect of variation in other parame- 
ters such as actuator or compensator dynamics or sensor location can also be studied. 
Primary subroutines called are PKEXT, MATIT, MNLOOP, DETIT, FINSCN, AT62, PREDICT, 
FUN, DEFINE, SCHEDUL, AERO, AEROSP, CXQR, MATINV, CXDEV, IUNI, and CSIUNI. Key 
parameters are defined as follows: 

t = 1 Use matrix iteration to obtain initial estimates of roots. 

= 0 Do not calculate initial roots using matrix iteration. 

I TRACE ( I ) Array of integers indicating whether to trace Ith root either from 

input or using IOPT1 or I0PT2 options. These parameters also 
indicate type of curve used in extrapolating for and s-j^. 

(See section "Determinant Iteration.") 

= 1 Linear. 

= 2 Qaadratic. 

= 3 Cubic. 
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KVAR 


IPKPLT 


NV 

IGAIN 

DV 

VO, HO, 
and RHOO 

DVEL + , DRHO + 
and DH + 

GN ( I , J ) 

DELGAIN(I, J) 

VV, H, 
and DENS 


Integer indicating type of variation to be performed. If negative, 
any corresponding plot output will be with respect to dynamic pres- 
sure. Can be reset interactively at plot time of IPKPLT > 0. 

= ±1 Velocity variation. 

= +2 Altitude variation. 

= ±3 Density variation. 

= 0 No plots. 

= +1 Root loci. 

= +2 G, oo plots, 
n 

= ±3 Both root loci and G / w plots. 

n 

= + Allows . alteration of titles, scales, curve deletion, etc., at 
plot time. 

Number of increments in parameter being varied. 

( = 0 No gain variation. 

= 1 Gain variation. 

Increment in parameter variation 

Initial values of velocity, altitude, and density 

r Incremental change in velocity, density, and altitude. Only one 

is set to DV, depending upon KVAR. Others are zero. 

Current value of gain on transfer function relating actuator I output 
to sensor J input. 

Increment in gain GN(I,J) for gain variation. 

Current values of velocity, altitude, or density during parameter 
variation. 


XVC and XV Current value of independent parameter being varied. 

A A 

SI I and S2I Array of initial estimates, s-j^ and s-j^ of roots to be traced. 

These need to correspond to initial value of parameter (velocity, 
altitude, density, or gain) variation. S2I, only, is used initially 
in MATIT if IOPT2 =1 as an estimate (not required) of roots to be 
traced. MATIT will return values for both SI I and S2I in this case. 
S2I is always assumed to be the better estimate in determinant 
iteration. 


t 


Internally 


generated parameter. 
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NF + 

ITEST + 


C0NFAC1 

C0NFAC2 


Number of characteristic roots for which to solve. 

Integer indicating whether there is a convergence problem in 
tracing a specified root. 

= 1 Iteration procedure converged. 

= 2 Iteration procedure failed to converge. 

Factor for converting input velocity to feet per second for use in 
subroutine AT6 2. 

Factor for converting input altitude to feet for use in subroutine 
AT6 2. 


C0NFAC3 Factor for converting density output by subroutine AT62 in slugs 

per cubic foot into units corresponding to input model. 

I PH and JPH Control (IPH) and sensor (JPH) pair for which phase error is to be 
introduced. 


PHASE 


Phase error, in degrees, to be introduced in IPH control, JPH 
sensor transfer function. 


IPRINT 


XMACH 

QBAR + 

KFIT 


Option indicating how often to print characteristic root data to 
screen (or to file OUTPUT) during interactive session; 

IPRINT = N, prints every Nth iteration. All data are written to 
file TAPE6, regardless of IPRINT. 

Mach number at which aerodynamic data were generated. 

Dynamic pressure corresponding to current velocity and density. 

Type of interpolation to be performed to obtain aerodynamic force 
at particular frequency when k aerodynamic tabular data are being 
used rather than p-plane aerodynamic approximations. 

If KFIT is negative, the interpolation will be for Re(Q(k)) 
and Q* (k ) (see unsteady aerodynamic force subsection) 

= ±1 Linear. 


= ±2 quadratic. 

= +3 Cubic spline. 


P1 + and P2 t Nondime nsionalized estimates, s± and si 2 of characteris tic 

root prior to determinant iteration at each value of independent 
variable. Upon returning from determinant iteration, P2 is 
current converged root, if found, and PI is preceding converged 
root. 


t 


Internally generated parameter. 
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P + Value of estimate for nondimensionalized root in MATIT corresponding 

to root in determinant iteration procedure for which convergence 
difficulty was incurred. Upon return from MATIT, P is replaced by 
actual root, if found. 


pii(d + 

and P22( I ) + 


P22(I) is value of converged Ith root X^ for current value of 
independent variable, and P1 1(1) is value of converged Ith root 
for preceding value of independent variable. P22(I) is employed 
to prevent iterations for Ith root from converging upon pre- 
viously determined root by defining determinant 

I 1 - 1 

DET = DET/ 11 (p - V) 


PP(I,IN) + Array of up to four preceding values of converged INth root. 

(1=1,4) and P11(IN) is equivalent to PP(4,IN). Array is used for predicting 

PPF(I) new root estimates PI and P2. 

NORD* Order of extrapolating curve currently being used to predict 

estimates PI and P2. 

STEPMX* Maximum step increase allowed in amplitude of characteristic root 

in determinant iteration. 


SIGN + 


NCHECK + 


ITRACK 1- 


NFINE 


IDEF 


t 


r < 0 Indicates change in stability has occurred in current root 
) being traced. 

V > 0 No change in stability has occurred. 


NF x 1 vector recording successive number of times noncovergence 
has occurred in roots being traced. If NCHECK ( IN ) > 5, the INth 
root is no longer traced. 


NF x i vector initially having 1 in all elements. If 

NCHECK ( IN) > 5, I TRACK (IN) is set to 0, and computations for INth 
root are bypassed 

Number of subintervals into which increment in independent 
variable is divided in performing finer scan for particular 
root. Number of values at which INth root is determined in 
finer scan. N = NFINE + 1. 

"= 1 Make precise determination of point where first change in 
stability occurs. 

.= 0 Bypass precise determination. 


t 


Internally generated parameter. 
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Vector array giving frequency, flutter speed, altitude, density, 
and damping obtained by interpolating for value of independent 
variable at which damping in INth mode is zero; i.e., at which 
change in stability occurred. 

IT* and NIT Number of iterations employed in seeking characteristic root, and 
maximum number allowed. 

ICUT* and NCUT If predicted step increases magnitude of determinant, step size 

is halved up to NCUT times in attempt to improve estimate of 
characteristic root. ICUT is number of times halving has 
occurred. 


EPSI If either relative or absolute change in estimate for characteristic 

root is less than EPSI, convergence is assumed. 

PROOT^ Estimate of characteristic root at which determinant is being 

eva luated. 


SCHEDUL 


I DAMP 


Complex number, for each control sensor pair, returned from 

SCHEDUL function subprogram giving contributions of scheduling 
component and any phase error. 



1=3 g. = s/o) g 

i / n . s . 
' 1 l 


r~ 0 Aerodynamic forces are computed as function of k = Im(p). 

= 1 Aerodynamic forces are obtained from p-plane approximation by 
using coefficients previously generated. 

ISPLANE / = 2 p-plane approximations to unselected set of modes is 

determined. 

= 3 p-plane approximations to selected subset of modes is 
^ determined. 


ICONS YS 


IDMULT 



Control system is included. 

No control system. 

Control system denominator is cleared. 

Control system denominator polynomial is not cleared. 


t 


Internally generated parameter. 
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= 1 Aerodynamic and inertial hinge moments are neglected in control 
rows of equations of motion 

ICSACT < 

= 0 Control surface is included as full degree of freedom, and 
effects of inertial and aerodynamic hinge moments are 
^ considered. 

Subroutine PKEXT . - In subroutine PKEXT , modal, sensor, aerodynamic, and control 
system data are read-in from data storage file TAP El ; the parameter variation loop is 
controlled; and characteristic root data are stored on data storage file TAP El for 
plotting later if IPKPLT t 0. Primary subroutines called are MATIT and MNLOOP. Key 
parameters used are I0PT2, ITRACE, IPKPLT, ISPLANE, ICONSYS, KVAR, IGAIN, NV, DV, 
DVEL, DH, DRHO, DELGAIN, SI I, S2I, NF, and IPRINT. 

Subroutines MATIT and MATITX . - Subroutines MATIT and MATITX, together, perform a 
matrix iteration to determine the characteristic roots. They are used to obtain 
initial estimates and °f t ^ e roots at the beginning value of the 

parameter (velocity, altitude, density, or gain) being varied if I0PT2 = 1 and/or 
to aid in obtaining a specific root when the determinant iteration process has con- 
vergence difficulties. (See appendix D. ) Primary subroutines called are AERO, 
AEROSP, SCHEDUL, CXQR, and MATINV. Key parameters used are VV, H, DENS, S2I, P, and 
ITEST. 


Subroutine MNLOOP. - In subroutine MNLOOP, the independent variable is incre- 
mented, initial estimates, of characteristic roots being traced are obtained, and 
(for each increment), the characteristic roots are determined using the determinant 
iteration method, if possible, or the matrix iteration method if the determinant 
iteration process fails to converge. Primary subroutines called are AT62, PREDICT, 
DETIT, MATIT, and FINSCN. Key parameters used are ITRACE, KVAR, IPKPLT, IGAIN, DVEL, 
DH, DRHO, DELGAIN, W, H, DENS, NF, ITEST, PI, P2, PI 1 , P22, PP, NORD, STEPMX, SIGN, 
NCHECK, ITRACK, NFINE, IDEF , IT, NIT, ICUT, NCUT, XMACH, C0NFAC1 , CONFAC2, C0NFAC3 , 
XVC, and IPRINT. 

Subroutine PETIT . - Subroutine DETIT performs the actual determinant iteration to 

obtain X , the INth characteristic root corresponding to P2* See section 
I N 

"Determinant iteration. " Subroutines called are FUN1 AND FUN. Key parameters used 
are ITEST, PI, P2, STEPMX, IT, NIT, ICUT, NCUT, EPSI, and PROOT. 

Subroutine FINSCN . - Subroutine FINSCN traces the IN characteristic root for 

smaller increments in the parameter being varied (velocity, altitude, density, or 

gain) in the same manner as MNLOOP in case of a change in stability or convergence 

difficulties. A new value is returned for P2 = X . The first time a change in 

IN 

stability occurs, the precise point where it occurs is obtained. Subroutines called 
are AT62, DETIT, MATIT, DEFINE, and PREDICT. Key parameters used are KVAR, IGAIN, 
DEVEL, DH, DRHO, DELGAIN, W, H, DENS, ITEST, PI, P2, PP, PPF, NORD, NFINE, IDEF, 
XMACH, C0NFAC1 , C0NFAC2, C0NFAC3, XV, and IPRINT. 

Subroutine AT62 .- Subroutine AT62 determines the density and speed of sound 
at a specified altitude. Key parameters used are H*C0NFAC2, DENS/CONFAC3, and 

v*confaci /xmach. 

Subroutine PREDICT . - Subroutine PREDICT fits a polynomial of order n < 3 to 
n + 1 points and evaluates the function at a specified value. It is used to extrap- 
olate for s-l and at the current value of the independent parameter (veloc- 

ity, altitude, density, or gain) based upon the values of the corresponding n + 1 
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converged roots. It is also used to interpolate for the precise point at which a 
change in stability occurs. Key parameters used are PI, P2, XVC, XV, PP, PPF, and 
NORD. 

Subroutines FUN1 and FUN .- In subroutine FUN1 and FIJN. the final matrix model is 
formed and the value of its determinant is obtained corresponding to the current 
estimate ^s INn+1 ^ of the root X^. Primary subroutines called are AERO, AEROSP, 

CXDEV, and SCHEDUL. Key parameters used are VV, DENS, PROOT, SCHEDUL, IDAMP, 

ISPLANE, ICONSYS, IDMULT, and ICSACT. 

Subroutine DEFINE .- Subroutine DEFINE uses PREDICT to interpolate for the pre- 
cise value of the independent variable (velocity, altitude, density, or gain) at 
which a change in stability occurs. (See section "Precise determination of point at 
which change in stability occurs.") Subroutine called is PREDICT. Key parameter 
used is V. 

Subroutine SCHEDUL .- Subroutine SCHEDUL provides for scheduled compensation of 
the signals from the sensors prior to sending them to the actuators. It also pro- 
vides for the inclusion of a phase error into the compensation for one control-sensor 
pair. See section "Compensation Cations" and appendix B. Key parameters used are 
PROOT, PHASE, IPH, JPH, XMACH, and QBAR. 

Subroutine AERO .- For ISPLANE = 0, aerodynamic forces as a function of 
k = Im(p) are determined in subroutine AERO. Subroutines called are IUNI and 
CSIUNI. Key parameters used are IM(PROOT) and KFIT. 

Subroutine AEROSP .- For ISPLANE ± 0, the aerodynamic forces are evaluated at 
p in subroutine AEROSP by using the p-plane approximations generated in overlay 
PPLANE. Key parameter used is PROOT. 

Subroutine CXQR .- The eigenvalues of a complex matrix are determined by subrou- 
tine CXQR along with selected eigenvectors. 

Subroutine MATINV .- Subroutine MATINV determines the inverse of a real matrix. 

Subroutine CXDEV. - Subroutine CXDEV computes the determinant of a complex 
matri x. 

Subroutine IUNI .- Subroutine IUNI performs linear or quadratic interpolation 
needed to obtain the aerodynamic forces at a particular value of k, depending upon 
KFIT = ±1 or ± 2 , respectively. 

Subroutine CSIUNI .- Subroutine CSIUNI performs cubic spline interpolation needed 
to obtain the aerodynamic forces at a particular value of p if KFIT = ±3. 


Overlay AEROPLT 

Plots are constructed which show how the oscillatory aerodynamic forces vary 
with reduced frequency and depict how the p-plane approximation fits the data. 
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Primary subroutines called are APLOT, MULTPT, IUNI, and CSIUNI. Key parameters are 
defined as follows: 


ISPLANE 


0 Plot of particular element of aerodynamic force matrix shows 
input data and curve defined by interpolating for points between 
these data. 

i 1 0 Additional curve showing the p-plane approximation is 
c displayed. 


MODE A ( I, J) 


'= 1 Plot of (I,J) aerodynamic force element is made. 
.= 0 (default) no plot is made for ( I , J ) th element. 


KMIN" 1 ’, KMAX* Minimum and maximum reduced frequencies for which aerodynamic 
forces are available. These set the range of reduced fre- 
quencies over which the data are plotted. 

IAPLT Indicates not only that aerodynamic data are to be plotted but 

also indicates which type of curves used to fit data are to be 
plotted. 


= 1 Plot interpolated curves, depending upon KFIT, and p-plane 
approximation if ISPLANE t 0. 


= -1 Plot only interpolated curves regardless of ISPLANE. 

KFIT Indicates type of interpolation to be performed in obtaining 

aerodynamic forces at a particular frequency. When KFIT is 
negative, kQ* (k) will be plotted for the imaginary part of 
the aerodynamic force (see the unsteady aerodynamic force 
subsection) . 


= ±1 Linear. 

= ±2 Quadratic. 

= ±3 Cubic spline. 


Subroutine APLOT . - KMIN and KM AX are determined in subroutine APLOT. Aero- 
dynamic forces at 50 intermediate values of k between KMIN and KM AX are also deter- 
mined by using interpolation and PPLANE approximating functions (if ISPLANE t 0 
and IAPLT > 0). Plot routine MULTPT is called. Subroutines called are IUNI, 

CSIUNI, and MULTPT. Key parameters used are all those defined for overlay AEROPLT • 


Subroutine MULTPT .- Subroutine MULTPT actually calls the plot labeling routines, 
line plot routines, axes routines, etc., and creates the plot vector file which 
defines the actual plots. These plot vector files can either be sent to the screen 
immediately for interactive plotting or simply saved and plotted after execution via 
any plotting device available. Subroutines called are library plotting routines for 
generating axes, drawing lines, labeling, etc. 


t 


Internally generated parameter. 
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Subroutine IUNI .- Subroutine IUNI performs linear or quadratic interpolation for 
the aerodynamic forces. 

Subroutine CSIUNI .- Subroutine CSIUNI performs cubic spline interpolation for 
the aerodynamic forces if KFIT = 3. 


Overlay PKPLOT 


Plots are constructed which portray the stability characteristics of the air- 
craft by using data generated in overlay PKFLUT if IPKPLT t 0 during its execu- 
tion. All available sets of data are listed if IPKPLT > 0. Also in overlay PKPLOT, 
the option of removing unwanted sets of characteristic root data is available. Pri- 
mary subroutines called are PKPLT, PLOT, RLPLT, SCALE, and DELPLOT. Key parameters 
are defined as follows: 

INTERAC* Indicates whether PKPLOT is being executed in interative mode. 

INTERAC is true if ICASE > 0. 

.T. All input is from file INPUT and output is to file OUTPUT. 

Ihese are usually assigned to the keyboard and screen, 
respective ly . 

.F. All input is from file TAPE2 , and output except plots is to 
file TAPE6. 


IPKPLT 


IDELPLT 


Integer indicating which type of characteristic root plots to 

generate and whether modifications to labels, scaling, and other 
plot parameters are to be allowed on a plot-per-plot basis. 

>0 Modifications are allowed. 

<0 Only the last set of data is plotted and no changes to plot 
parameters are allowed during generation of plots. 

= +1 Generate plots of root loci versus velocity, altitude, or 
density as determined by KVAR =1, 2, or 3, respectively, 

or versus corresponding dynamic pressure if KVAR < 0, or 
versus gain if IGAIN = 1 during root determination. 

= ±2 Generate plots of C and U) versus velocity, altitude, or 
density depending upon KVAR during root determination (or corre- 
sponding to dynamic pressure if KVAR < 0) or versus gain if 
IGAIN = 1 during root determination. 

= ±3 Generate both ) and root loci plots. 

n 

= 1 Delete unwanted characteristic root data from data storage 
file TAPE1 . 

= 0 Do not delete any data. 


t 


Internally generated parameter. 
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I COUNT 

IPLT and 
IPLT2 

I COMBI N 


PTITLE 


IDASH 


I SAME 


KVAR 


Number of sets of characteristic root data to be plotted concurrently 
out of total number of sets of data generated during different 
cases. 

Array of numbers indicating which sets of characteristic root data 
are to be plotted or saved, respectively. 

Integer indicating whether to combine plots of different sets of 
characteristic root data into one plot (all on same page) or to 
generate plots for separate pages. 

= ^ Combine plots. 

= 0 Do not combine plots. 

80-character description of data being plotted. When plots are 
being combined, this will appear on plot being generated in lieu 
of the original plot titles. 

Indicates type of line used to connect points in data currently 
being plotted. This may be changed for each set of data being 
plotted. 

= 0 (solid line) 

= 1 _ __________ 

= 2 


= 3 


= 4 


= 5 


= 6 


= 7 


Integer indicating how many of next sets of characteristic 
root data are to be plotted by using same PLOTP namelist 
data as currently being read-in, including current set. 
(Default = 1). Program will not request another PLOTP 
namelist to be input until n = ISAME sets of data have been 
plotted. 

Sign of KVAR indicates whether data are to be plotted versus 
velocity, altitude, or density (KVAR > 0) or versus dynamic 
pressure (KVAR < 0) . 
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MODEPK(I) Integer indicating whether to generate a plot showing Ith character- 

istic root variation, 

= 1 Generate plot of data (Default) 

= 0 Bypass Ith root information. 

XRLORG, XRLEND, Lower limit, upper limit, increment 

XRLDEL, and number of minor divisions between 

NDXRL in root-locus plots. (Optional) 

YRLORG, YRLEND, Lower limit, upper limit, increment between major tic marks, and 
YRLDEL, and number of minor divisions between major tic marks for ordinate 

NDYRL in root-locus plots. (Optional) 

VORG, VEND, Lower limit, upper limit, increment between major tic marks, and 

VDEL, and NDV number of minor divisions between major tic marks for abscissa 

in C and co plots. (Optional) 
n 

FORG, FEND, Lower limit, upper limit, increment between major tic marks, and 

FDEL, and NDF number of minor divisions between major tic marks for ordinate 

in co plots. (Optional) 
n 

GORG, GEND, Lower limit, upper limit, increment between major tic marks, and 

GDEL, and NDG number of minor divisions between major tic marks for ordinate 

in C plots. 

Integer indicating whether to continue plotting. This option 
allows user to generate additional plots with characteristic 
root data that either were not included before or for which 
changes in plots are desired. 

= 1 Continue plotting. 

= 0 Return to main program and continue with new case or terminate 
execution. 

NFR^ Total number of roots traced in current data set being plotted. 

NVR f Total number of increments in independent variable for current 

data set 

IDYNAM^ Indicates whether to convert to dynamic pressure rather than plot 

original independent variable. 

= 1 Convert to dynamic pressure. 

= 0 Do not convert. 


internally generated parameter. 


I OP 




between major tic marks, and 
major tic marks for abscissa 
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Subroutine PKPLT. - Subroutine PKPLT actually reads in the plot data, converts 
the independent variable to dynamic pressure, if desired, and calls plotting routines 
to generate desired plots. Subroutines called are PLOT and RLPLT. Key parameters 
used are IPKPLT, NFR, NVR, MODEPK, and IDYNAM. 

Subroutine PLOT - Entry points PL0T1 and PL0T2 .- Subroutine PLOT generates the 
plots of £ and oo n versus the corresponding independent variable desired (veloc- 
ity, altitude, density, dynamic pressure, or gain). Subroutines called are SCALE and 
library plotting routines for generating axes, drawing lines, labeling, etc. Key 
parameters used are ICOMBIN, IDASH, MODEPK, NFR, NVR, VORD, VEND, VDEL, NDV, FORG, 
FEND, FDEL, NDF, GORG, GEND, GDEL, and NDG. 

Subroutine RLPLT . - Subroutine RLPLT generates the root locus plot. Initial 

points are identified with the mode number. Every 10th point is identified with a + 

for nongain variation plots. Gain points are identified on the gain variation plots, 

which are 0, .5, 1., and 2. times the nominal gain incorporated by input into 

(T*) . Subroutines called are IUNI, SCALE, and library plotting routines for 

L i,j 

generating axes, drawing lines, labeling, etc. Key parameters used are ICOMBIN, 
IDASH, MODEPK, NFR, NVR, XRLORG , XRLEND, XRLDEL, NDXRL, YRLDRG, YRLEND, YRLDEL, and 
NDYRL. 


Subroutine SCALE . - Subroutine SCALE determines an appropriate power of 1 0 to be 
used for tic marks on axes which are 2., 4., 5., or 10. in. in length and returns the 
corresponding scale factor required to scale the data. 

Subroutine DELPLOT. - Subroutine DELPLOT deletes unwanted sets of characteristic 
root data, generated by using overlay PKFLUT, from data storage file, TAPE1 . Key 
parameter used is IPLT2. 
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INPUTS FOR SAMPLE CASES 

The TAPE 2 inputs required to generate the plots shown in figures 8 to 1 6 are 
presented in this appendix. These inputs exercise the major options of STABCAR. 
They are presented to serve as a guide for users to follow in preparing input data. 
Reference should be made to table II for the definitions of the input variables. 


Cases 1 to 4 (Figs. 8 to 1 1 ) 

The inputs for sample cases 1 to 4 (figs. 8 to 1 1 ) are presented as a unit in 
table GI. These cases are run in a batch mode (ICASE < 0). Case 1, which generates 
the data for figure 8, has ICASE = -2. This means that TAPES (table III) must be 
local in order to make the airplane generalized aerodynamic forces available. In 
cases 2 to 4, ICASE = -1 means that a batch mode of operation is in effect and all 
data must be obtained from TAP El (table V) , which was generated in the previous case 
and from TAPE2. Note that in case 3 the array S2I (s-j^) has been set to zero. This 

was done because, with IPS = 0 (the default), S2I contains the roots found at the 
end of the last case having IFLUT = 1. That set of roots is for h = 0 (high 
dynamic pressure), whereas case 3 begins at a low dynamic pressure. Setting S2I 
to zero restarts the estimation process; thereby, convergence difficulties are 
avoided that would have occurred since S2I is used in starting the matrix iteration 
(I0PT2 = 1) search for initial estimates. (See appendix D. ) These cases could have 
been run interactively with the plots being returned during job execution by making 
ICASE > 0 and loading LIBFTEK and Tbchtronix interface plotting routines, which 
support interactive graphics, or their equivalent. 


Case 5 (Fig. 12) 

Tables GII, GUI, and GIV contain the inputs required to generate figure 12. 

The modal deflections at two accelerometer locations are needed to implement 
the flutter control law. TAP El 0 (table IV) which contains the mode shape data, 

TAPES, and TAPE2 are made local. For clarity a new TAPE2 is generated? it is shown 
in table GII. Alternatively, one could have modified the initial TAPE2 corresponding 
to case 1 . 

The modal deflections determined in this run are differenced and this difference 
is read onto TAPE5. None of these operations are shown. Next a run is made which 
defines the control law. The inputs are shown in table GUI. 

At this point the user saves TAPE1 which now contains the control-law dynamics 
and selected sensor data. This TAPE1 will be employed not only to complete case 5 
but also as part of the input for cases 6, 8, and 9. Thble GIV presents TAPE2 input 
required to complete case 5. These runs could have been done interacti ve ly by set- 
ting ICASE = 1. Note that ISCHEDUL (3,1) = 1 which means that the scheduled con- 
trol law factor defined in appendix B and figure 7 will be employed. 
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Case 6 (Fig. 13) 

The input required to generate figure 13 is now shown. The TAPE1 employed is 
as it existed prior to making the last two runs in case 5. Thus, the aircraft is 
modeled including two rigid-body, six elastic, and one control mode, the correspond- 
ing sensor deflections, and the definition of the flutter suppression system control 
law. ICASE is set to +1 in any convenient TAPE2. For ICASE = +1 none of the rest 
of TAPE2 is employed. The values of all parameters are contained on the TAPE1 which 
must be local. The interactive session then proceeds as shown in table GV. 


Case 7 (Fig. 14) 

The interative session during which figure 14 was generated is shown in 
table GVI. This figure shows the effect of phase errors in the control system upon 
the elastic mode gain root loci. The TAPE1 created in case 6 is employed as the 
primary source of input. Three runs are made. The first two of these runs deter- 
mine stability characteristics by using the p-p method after phase errors have been 
introduced. The third run combines the root loci plots from these two runs with 
phase errors with the corresponding plot from case 6 which had no phase errors. 


Case 8 (Fig. 15) 

Figure 15 presents gain root loci for the control-law poles (see fig. 7) and one 
of the open-loop actuator poles (s = -401 rad/sec). The loci are shown at the open- 
loop flutter point for N Ma = 0.86. Thble GVII shows the interactive session which 
created the figure. The TAPE1 saved prior to running the last two segments of case 5 
is the primary source of input. Modifications are made with the TAPE2 shown at the 
top of the table. 

The initial segment of this session traces one elastic root, mode 8, one actua- 
tor root, mode 13, and three control-law roots from zero gain to a gain of 0.07. The 
scheduled pole is not traced because this denominator is not cleared and, for a gain 
of zero, would result in division by zero in equation (2). Note that IDMULT = 1 in 
TAPE2 which means that the nonscheduled control-law denominator, actuator denom- 
inator, and sensor denominators are cleared in equation (2). The card image contain- 
ing the three zeros tells the program, when IPKPLT = 1, that no plots are to be 
generated in this initial batch (ICASE = -1) run. This first segment is generated 
knowing in advance that a series of runs will be required with the resulting plots 
combined into one composite plot. The mode number identifiers (fig. 15) are only 
placed on the first plot set. Consequently, one plot set - they can be plotted in 
any order - should contain all the modes starting from the initial value of the inde- 
pendent variable. The +1 ending TAPE2 makes ICASE > 0 so that the interactive mode 
of operation is invoked. 

A series of interactive runs are made to generate the gain loci by using the 
p-k method. Note the first interactive run where IPS = 1. This value for IPS means 
that the initial estimates at the start of the last execution with IFLUT = 1 will 
be employed as initial estimates in the current run. Note also the runs where 
root 1 2 is being determined as a function of gain. This is an illustration of the 
case where two real roots combine to form a complex pair. The next to last run is a 
batch run where roots 11 and 14 are found with the p-p method. Finally, the last 
case combines the plot sets to form figure 1 5. 
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Case 9 (Fig. 16) 

The inputs required to generate figure 16 are shown in table GVIII. Figure 16 
is an open-loop locus of the short period root with altitude. Three interactive runs 
are made. The TAP El generated in case 2 is made local and ICASE is set to +1 in a 
TAPE2 • The first two runs obtain the root locus for the p-k and the p-p method. The 
third run combines the two plot sets to produce figure 16. 
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TABLE GI.- INPUT FOR SAMPLE CASES 1 TO 4 


-2 

CASE 1. P-PLANE FIT AND AERODYNAMIC FORCES PLOTS. (FIG 8) 
$ INPUT NM= 1 6 , NR== 2 , NC= 4 , ICSACT= 1 , NK=1 2 ,KFIT=2 , 


ISPLANE=3 , NCOEF=7 , NKK( 1 )=2*5 , IC0F( 1 , 1 )=0 ,0 , 1 , 0 , 1 , TCOF( 1 , 2)= 1 , 1 , 
NCOL(1)=0,BN(1,1)=0. 0939 ,0.188, 0.281 6, 0.376, 

CSURFEF( 3)= 1 . , I SELECT* 1 , IAPLT= 1 ,C= .5961 38, 

M0DEA( 3 , 3 ) = 1 , MODEA( 5 , 3 ) = 1 , MODEA( 8 , 3 ) = I , 

MODEA( 3 , 5)= 1 ,M0DEA( 5 , 5) = 1 ,MODEA( 8 , 5)= 1 , 

MODEA( 3 , 8 ) = 1 , M0DEA( 5 , 8 ) = 1 , M0DEA( 8 , 8 ) = 1 , 

MODEA( 3 , 15)= 1 ,MODEA( 5 , 1 5) = 1 ,M0DEA( 8 , 1 5)= 1$ 

$SELECT NMODE S= 9 , NRNEW= 2 , NCNEW= 1 ,N0C( 1 )* 1 , 2 , 3 , 4 , 5 , 6 , 7 ,8 , 1 5 $ 

-1 

CASE 2. OPEN LOOP STABILITY CHARACTERISTICS VERSUS ALTITUDE ,DV=0 . 5 (FTG 9) 
$ INPUT T.APLT=0 , I SPLANE=0 , IFLUT* 1 , NV=60 , ICSACT* 1 , 

T.SELECT=-1 , 


DV=0 • 5 
XMACH=0 .86 
NIT=26 
I0PT2=1 
VORG=0. 

FORG=0 
G0RG=-1. 
XRLORG=-300. 
YRLORG-O. 
I.TRACE( 1 ) = 8* 1 


, H0= 30. 

,C0NFAC1=3. 28084 
,NCUT>25 
,IPRINT=5 
,VEND=30 • 

, FEND* 60 
, GE ND= 1 . 0 
,XRLEND=100. 
,YRLEND=400. 

$ 


,KVAR=2 

,CONFAC2=3280.84 

,NFINE*4 

,IPKPLT=-3 

,VDEL*10 

,FDEL=20 

,GDEL=0.5 

,XRLDEL=100. 

, YRLDEL* 100. 


,EPSI= 1 .E-6 , 
,CONFAC3=515. 3787 
, I DAMP* 2 
,MODEPK( 1 )=8* 1 
,NDV= 1 0 
,NDF=4 
,NDG=5 
,NDXRL=5 
,NDYRL=4 


$SELECT NMODES=8,NRNEW=2,NCNEW=0,NOC(1)=1,2 > 3 > 4 > 5,6,7 ,8$ 

-1 

CASE 3. OPEN LOOP CHARACTERISTIC ROOTS VERSUS DENSITY, DV=0 .025 (FIG 10) 
$INPUT NV*28 , DV-0 .025 ,RHO0=0.05 ,V0=269.64 , S2I=25*( 0 , 0) , 

KVAR* 3 ,VORG=0 .0 ,VEND=0.80 ,VDEL=0.20 , 

NDV=10 ,GORG=-0 • 2 ,GEND*0.6 ,GDEL=0.2 , 

NDG* 10 ,XRL0RG=-1 20 ,XRLEND=20,XRLDEL=20 , 

NDXRL=4 $ 


-1 


CASE 4. OPEN LOOP STABILITY CHARACTERISTICS VERSUS DYNAMIC PRESSURE. (FIG 11) 
S INPUT IF L13 T=0 , KV AR=- 3 ,VORC=0. , VEND=30000 . , VDF,L=5000., NDV*10 $ 

0 


TABLE GII MODAL DEFLECTIONS AT ACCELEROMETER LOCATIONS 


-2 

CASE 5A. MODAL DEFLECTIONS AT ACCELEROMETER LOCATIONS 
$INPUT NM=16,NR=2,NC=4,NS=2,IFLUT=0,IPKPLT=0, 

ISENSEf 1 , TNTERP* 1 , ISYM* 1 , NPMX=60 , 1 PLATE( 1 )= 1 , 

NSECTNS=3 ,NNODES=100,ISS(1,1)=1 ,60,61 ,92,93,100, 

XCG=6. 71576,X0(1)=6. 4470,0.0,9. 1 328, YO (3) =. 25 15 , 
R0(1)=1. 7799, 1.570796,. 5 175 $ 

$SENLOC XS(1)=7.1438,7.3363,YS(1)=2.083,2.134, 
ITYPE(1)=1,1,NSS(1)*1,1 $ 

0 
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TABLE GUI.- SELECTION OF SENSOR DATA AND DEFINITION OF CONTROL LAW 


-2 

CASE 58. SELECTION OF SENSOR DATA, AND DEFINITION OF CONTROL LAW. 
$INPUT NM= 16, NR= 2 , NC= 4 , TC SAC T= 1 , NK= 12 1 KFIT=2,NS=l, 

I SPLANE= 3 , NCOEF= 7 , NKK( i ) = 2*5 , ICOF( 1 , 1 )=0 ,0 , 1 ,0,1 ,IC0F( 1 ,2) = 1 ,1 , 
NCOL( 1)=0,BN( 1, 1)=0. 0939, 0.188, 0.2816, 0.376, 

CSURFEF( 3)= 1 • , I SELECT= 1 , TAPLT=0 , C= . 596 1 38 , 

MODEA( 3 , 3)= 1 ,MODEA( 5 , 3)= 1 ,MODEA( 8 , 3)= l , 

MODEA( 3 , 5 ) = 1 , MODEA( 5 , 5) = 1 , M0DEA( 8 , 5) = 1 , 

MODEA( 3 , 8)= 1 , MODEA( 5 , 8 ) = I , MODEA( 8 , 8)= 1 , 

MODEA( 3, 15)= 1 ,MODEA(5, 15)= 1 ,MODEA( 8 , 1 5) = 1 , 


, H0=4 .572 
,CONFACl=3. 28084 
, NCUT= 25 
, IPRI NT=5 
, VEND=30. 

, FEND=60 
,GEND= 1.0 
,XRLEND=100. 

, YRLEND-400. 


,KVAR=2 

,C0NFAC2= 3 280. 84 
, NFTNE=4 
, I PKPLT=-3 
,VDEL=10 
,FDEL= 20 
,GDEL=0 • 5 
,XRLDEL“ 1 00. 
,YRLDEL= 100. 


, EPSI= 1 . E-6 , 

, CONE AC3=515. 3787 

, I DAMP =2 

,MODEPK( 1 ) = 8* 1 

,NDV=10 

,NDF=4 

,NDG=5 

,NDXRL= 5 

,NDYRL=4 , 


DV-0.5 ,NV= 1 
XM AC H= 0 • 8 6 
NIT=26 
IOPT2= 1 
VORG=0 . 

FORG=0 
GORC=- 1. 

XRLORG=-*300 . 

YRL0RG=0. 

TTRACE( 1 ) = 8* 1 , 

I C0NSYS= 1 ,GN(3,1)=0. , ISCHDUL( 3 , l )= 1 , 

I SENSE=0 , IFLUT=0 ,IPKPLT=0 ,INTERP=0 , ISYM=1 ,NPMX=60 ,TCSREAD=1, 

NSECTNS= 3 , NN0DES= 1 00 , I S S( 1 ,1)=1 ,60,61,92,93,100, TPLATE( 1 ) = 1 , 

XCG=6. 7 1576, X0(1) =6. 4470,0.0,9. 1328, 

Y0(3)=0. 2515, R0 ( I ) = 1 . 7799, 1 . 570796 , . 5 1 75 $ 

$SENLOC XS( 1) = 7.1438, 7. 3363, YS( 0 = 2.08 3,2.134, 

I TYPE( 1) = 1»1 , NSS( 1 )= 1 , I $ 

$SENINP I SDYN=0 ,I0RD(1)=2 S 


3 1 

$F T LTTN IF I LTER=6 ,6 , 6 , 6 , 1 , 

AFN( 1 , 0 = 250. , AFD( 1 ,0 = 250. , 1 . , 

AFN(1 ,2)=2.079E5,220. ,1. , AFD( 1 , 2)=2 . 704E5 , 400 . ,1 . , 

AF N ( 1 , 3 ) = 3 . 8 4 1 6 E4 , 4 0 . ,1. , AFD( 1 , 3) = 3.86 1 2E4 , 100. ,1. , 

AFN( 2 , 4 ) = — 1 .1 554 , AFD( 1 ,4) = 2. ,1. $ 

0 0 

$ACTINP IACT( 3) = 1 , AACTN( 1 , 3) = 6 . 583E 1 4 , 

AACTD( 1 ,3)=6.583E14,2.348E12,2.9939E9,3.607E6,1.7373E3, 1. $ 

$ S ELECT NMODES=9 , NRNEW=2 , NCNEW= 1 , NOC( 0=1,2,3,4,5,6,7,8,15$ 

0 


TABLE GIV.- INPUT FOR SAMPLE CASE 5 


-1 

CASE 5C • OPEN LOOP STABILITY VERSUS VELOCITY (REDUCED ORDER MODEL, FIGURE 12A) 
$ INPUT I SELECT=- 1 , ICSREAD=0 ,KVAR= 1 , IFLUT= l , I SPLANE=0 , 

RHO0=0.6 101 ,V0=50. ,NV=75 , DV=4 , lAPLI^O , TPS=0 , S2T=2 5*(0 ,0) , 

TPKPLT=-3 , VORG=0 • , VEND=350 , VI)EL=50 • , 

NDV=10 ,GORG=-l . ,GEND= 1 .0 ,GDEL=0. 5 , NDG=5 $ 

$ SELECT NM0DES=4 ,NRNEW=0 , NCNEW= 1 , NOC( 1) = 3 , 5 , 8 , 1 5 $ 

-1 

CLOSED LOOP STABILITY VS VELOCITY (REDUCED ORDER MODEL, FIGURE 12B) 

$ INPUT GN( 3 , 1 )= 1 • , V0=50. , I PS=0 , S2I=25*(0 ,0) $ 

0 
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FIELD LENGTH OVERLAY (5.0) IS 100464 F'K FLUTTER PLOTS 

■ DO YOU UISH TO CONTINUE PLOTTING? (1/0)* 


APPENDIX G 


TABLE GVI II INTERACTIVE EXECUTION OF SAMPLE CASE 9 


TYPE IN TITLE 

? CASE 9. LOCUS OF SHORT PERIOD ROOT WITH ALTITUDE <P-X, FIGURE 16) 

TYPE IN /INPUT/ 

? $ INPUT I S P L ANE = 0 , IFLUT= 1 , IPKPLT=- 1 , XV AR = 2 , H 0 = 1 6 , DV= 1 , NV =10, I P R I NT = 1 0 , 

? I TRAC E = 0 ,1,6*0, 10 PT 2=1$ 

1 DETERMINATION OF STABILITY CHARACTERISTICS - - > ST AB C AR < - - DATE: 83/09/12 

TIME : 11.33.04. 

CASE 9 LOCUS OF SHORT PERIOD ROOT WITH ALTITUDE (P-K, FIGURE 16) 

MACH #= .8600 


FIELD LENGTH-OVERLAY ( 4 , 0 ) IS 104264 P-K FLUTTER ANALYSIS 

INITIAL ROOTS DETERMINED BY MATRIX ITERATION 


FIELD LENGTH INCREASED TO 104320 
1THIS CASE IS FOR 


AN 

INITI AL 

DENSITY 

OF 

. 1664707E+00 


AN 

INITIAL 

VELOCITY 

OF 

. 2 5 3 759 7E + 03 


AN 

INITIAL 

ALTITUDE 

OF 

. 1 6 0 0 0 0 0 E+ 02 , 

AND 

AN 

INITI AL 

DYNAM I C 

PRESSURE OF 5359857E+04 


DO NOT TRY TO DETERMINE ROOTS CORRESPONDING TO GAIN SCHEDULING DENOMINATORS 

FOR G A I NS = 0 . 

0 INITIAL SI S2 

MODE 

2 ( - . 33 8 06 44 3E+00 , . 2 9 5 2 5 5 5 1 E + 0 1 ) ( - . 3 4 1 4 7 9 2 2 E + 0 0 , .298 2 

3 7 8 9E + 0 1 > 

# - ITER ROOT NBR ALTITUDE ZETA OMEGA-N 

ROOT PRESSURE 

2 2 16000G0E+Q2 1137569E+00 4777617E+00 ( 

- . 3 4 1 4 8 3 0E + 0 0 , . 2 9 8 2 3 7 9 E + 0 1 ) 5359857E+04 

4 2 . 6000000E+01 2 0 7 8 0 2 9 E 4- 0 0 .9194466E+00 ( 

- 1 2 0 0 4 8 9E + 0 1 , 5 6 5 0 9 4 4E + 0 1 ) 2444549E+05 

CURRENT PLOT TITLE IS: 

CASE 9 LOCUS OF SHORT PERIOD ROOT WITH ALTITUDE (P-K, FIGURE 16) 

HIT RETURN OR TYPE' SAME , IF OK; OTHERWISE TYPE IN NEW TITLE. 

? SAME 

DO YOU WISH TO GENERATE PLOTS OF STABILITY CHARACTERISTICS AT THIS TIME? 

? 0 

CONTINUE 
? 1 

TYPE IN TITLE 

? CASE 9. LOCUS OF SHORT PERIOD ROOT WITH ALTITUDE (P-P, FIGURE 16) 

TYPE IN /INPUT/ 

? 5INPUT ISFLANE=1 , IPS = 0 , SI I = 25*(0 , 0) ,52 I = 2 5 * ( 0 , 0) , IOPT2=l$ 

1 DETERMINATION OF STABILITY CHARACTERISTICS - - > ST A B C AR < - - DATE: 83/09/12. 

TIME : 1144 53 . 

CASE 9. LOCUS OF SHORT PERIOD ROOT WITH ALTITUDE (P-P, FIGURE 16) 

MACH t= 8600 


FIELD LENGTH-OVERLAY ( 4 , 0 ) IS 101750 P-K FLUTTER ANALYSIS 

INITIAL ROOTS DETERMINED BY MATRIX ITERATION 

UNSTABLE ROOT 


FIELD LENGTH INCREASED TO 102004 
1THIS CASE IS FOR 

AN INITIAL DENSITY OF 1664707E+00 

AN INITIAL VELOCITY OF 2537597E+03 

AN INITIAL ALTITUDE OF 1600000E+Q2, AND 

AN INITIAL DYNAMIC PRESSURE OF 5359857E+04 
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APPENDIX G 


TABLE GVIII.- Concluded 


DO NOT TRY TO DETERMINE ROOTS CORRESPONDING TO GAIN SCHEDULING DENOMINATORS 

FOR G A I NS = 0 

0 INITIAL SI S 2 


MODE 

2 ( - 4 1 972234E+00 , . 2 7 7 4 1 7 5 0 E + 0 1 ) ( -42396196E+00, 2 8 0 2 


1 9 7 0 E + 0 1 ) 

# - ITER ROOT NBR ALTITUDE 

ROOT 

1 2 1600000E+02 

- 4 2 3 96 1 9E+ 00 , 2802197E+01) 

3 2 6000000E+01 

- 1 6 1 3 7 8 3 E + 0 1 , 5159008E+01 ) 

CURRENT PLOT TITLE IS 


ZETA 

PRESSURE 

1 495938E+00 
5359857E+04 

2985435E+00 

2444549E+05 


CASE 9 


LOCUS OF SHORT PERIOD ROOT WITH ALTITUDE (P-P, FIGURE 


OMEGA-N 
45 10590E+Q0 
8603155E+00 


1 6 ) 


< 

< 


HIT RETURN OR TYPE SAME, IF OK; OTHERWISE TYPE IN NEW TITLE 
’ SAME 


DO YOU WISH TO GENERATE PLOTS OF STABILITY CHARACTERISTICS AT THIS TIME’ 
’ 0 


CONTINUE 
’ 1 

TYPE IN TITLE 

’ CASE9 LOCUS OF SHORT PERIOD ROOT WITH ALTITUDE (FIGURE 16) 

TYPE IN /INPUT/ 

’ $INPUT IFLUT=0, IPKPLT=1$ 

1 DETERMINATION OF STABILITY CHARACTERISTICS - - > ST A B C AR < - - DATE; 83/09/12 

TIME 11 50 35 

CASE? LOCUS OF SHORT PERIOD ROOT WITH ALTITUDE (FIGURE 16) 

MACH 8600 

DO YOU WISH TO GENERATE PLOTS OF STABILITY CHARACTERISTICS AT THIS TIME’ 

’ 1 

PLOT SET TYPE OF VARIATION # MODES PLOT TITLE 


1 ALTITUDE 0 OPEN LOOP STABILITY CHARACTERISTICS VER 

SUS ALTITUDE, DV = 0 5 (FIG 9) 

2 ALTITUDE 1 CASE 9 LOCUS OF SHORT PERIOD ROOT WITH 

ALTITUDE (P-K, FIGURE 16) 

3 ALTITUDE 1 CASE 9 LOCUS OF SHORT PERIOD ROOT WITH 

ALTITUDE (P-P, FIGURE 16) 

H TYPE IN NUMBER OF PLOT SETS TO BE PLOTTED," 

" CORRESPONDING PLOT SET #'S, AND WHETHER TO COMBINE THEM' 1 
’ 2 
2 3 1 

CURRENT PLOT TITLE IS: 

CASE9 LOCUS OF SHORT PERIOD ROOT WITH ALTITUDE (FIGURE 16) 

HIT RETURN OR TYPE SAME, IF OK; OTHERWISE TYPE IN NEW TITLE. 

’ SAME 

TYPE IN / PLOTP/ -KVAR , I DASH , MODEPK , ETC FOR PLOT SET # 2 

’ $ PLOTP I D ASH = 0 , XRLORG=-2 , X R L END = 0 , X RL D E L = 1 ,NDXRL = 2, 

’ YRLORG= 0 , YR L END = 6 , YRLDEL = 2 ,NDYRL = 4$ 


FIELD LENGTH OVERLAYS, 0> IS 0 7 7 2 5 6 PK FLUTTER PLOTS 

TYPE IN /PLOTP/ -KVAR , IDASH , MODEPK , ETC FOR PLOT SET # 3 

’ 5 P L OTP I DASH= 1 $ 


FIELD LENGTH OVERLAY(5,0) IS 077256 PK FLUTTER PLOTS 

” DO YOU WISH TO CONTINUE PLOTTING’ (1/0)" 

’ 0 

CONTINUE 
’ 0 
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TABLE X.- OUTLINE OF DATA ON TAPE2 AND/OR TAPE3 (INPUT) 




(a) 


(a) 


(a) 


ISENSE = 1 


Variable or 
namelist name 


ICASE Free-field integer 


HEADER Eight-word (A10) alphanumeric array 


$ INPUT. . .$ I Namelist 


$SENLOC • • • $ Namelist 





5 ICONSYS * 1 & ICASE = ±2 l $SENINP 
or 

ICHSEN = 1 & ICASE = ±1 


6 ICONSYS = 1 & ICASE = ±2 
or 

ICHFIL = 1 & ICASE = +1 


7 ICONSYS = 1 & ICASE = ±2 $FILTIN...$ Namelist 
or 

ICHFIL = 1 & ICASE = +1 


Repeat sets 6 and 7 for each control-sensor pair; terminate with 0,0 


8 ICONSYS = 1 & ICASE * ±2 $ACTINP...$ Namelist 
or 

ICHACT = 1 & ICASE = ±1 




I SELECT * 0 


IDELPLT * 0 


IPKPLT > 0 


IPKPLT > 0 


S SELECT 


IPLT2 


ICOUNT 


IF 


Namelist 


I Free -fie Id integer array 


Free-field integer 


Free-field integer array 
(Size = ICOUNT) 


14 IPKPLT > 0 & ICOMBIN = 1 | PTITLE 


ICOMBIN 

Free-field integer 

PTITLE 

Eight-word (A10) alphanumeric array 


IPKPLT < 0 


$PL0TP. ..$ 



Repeat set 1 5 up to ICOUNT times (see definition of ISAME) 


I OP Free-field integer 


Repeat sets 10 to 16 as desired, terminate with IOP = 0 
Repeat sets 1 to 1 6 as desired for each case, end with ICASE 
Required in all cases. 























































TABLE II.- DESCRIPTION OF TAPE2 INPUT DATS 


When required 


Algebraic 

variable 


FORTRAN variable 


Description 


Integer indicating type of case to be run 
{initialization or continuation) and which 
version of proqram (interactive or noninter- 
active) to invoke. 

Initialization case. All data are (re)initial- 
ized by namelist and option parameters from 
TAPE2 and all array data (such as aerody- 
namics, generalized masses, and mode shapes) 
from TAPES and TAPE1 0. 

Continuation case. Case being run will use 

data previously stored on TAPE1 , data storage 
file created and used by program during 
execution. Changes only to namelist and 
option parameters are read-in from keyboard 
if I CASE = + 1 and from TAPE2 if -1. Changes 
to array data are made according to I SELECT. 

Execution stopped. 

Alphanumeric description of case being run. 

Total number of modes: rigid, elastic, and 


Number of rigid-body modes (see section on 
"Sensor input definitions"). 

Number of control modes. 

Number of reduced frequencies at which aerody- 
namic data have been computed. 

Number of gust forces included as columns in 
aerodynamic array input. Provision has beer 
made to read-in aerodynamic qust forces whic 
may be included as additional (NG) columns i 
aerodynamic data. These data are not used t 
the proqram. (Default = 0) 

Number of sensors. (Default = 0) 


Rewind TAPES input file at beginning of case. 
This a II ows for reinitialization of original 
array data. (Default.) 

Do not rewind TAPES. This allows new aerody- 
nami c and structural data to be read-in. 


Rewind TAP K2 input file immediately after rent 
ing namelist INPUT, allowing rei ni t.i a 1 i zati c 
of sensor and control system data. 

Do not rewind. (Default) 

Indicates whether model definition arrays are 
to be altered for current case. Different 
modal combinations can be selected, differ**] 
subsets of aerodynamic data in terms of Fre 
quenuy can be selected, etc. Refer to def i 
nition of namelist. SELECT. 

Model definition arrays are read from TAPES a 
data are selected from this input. Refer ti 
definition of IKEWb if reinitializing. 


Individual parameters are optional if default values are accept able 



ICASE = +2 


CSURFEF ( I ) 


2b or c 

C 


n u 

Ma 

XMACH 


Model definition arrays are read-in from TAPE1 
(a data storaqe file created and used by pro- 
gram during execution). These data represent 
model as it was at end of preceding case, 
incorporating any modal and/or frequency 
selections up to that point. New modal and/ 
or frequency selections are from this input. 


Model definition arrays will not be altered for 
this case. (Default) If ICASE = ±1, the 
arrays will be precisely those in effect at 
the end of the preceding case. 


Multiplier which modifies column of aerodynamic 
force matrices corresponding to ith control 
mode. The multiplication is only performed 
when ICASE - ±2 or I SELECT > 0. 

(Default =1.0) 


Reference chord length 


Mach number at which aerodynamic forces were 
generate d. 


+1 


- ±2 
- +3 
I STIFF 


# 0 
= 1 


Type of fit to be employed when interpolating 
input aerodynamic force data to obtain aero- 
dynamic forces at a particular frequency. " 
Used when ISPLANE = 0. When K FIT < 0, 
the interpolation is for Re(Q(k)) and 
O' (k) (see section "Unsteady aerodynamic 
forces " ) . 


Linear fit. 

Quadratic fit. (Default) 
Cubic fit. 


Integer indicating hew stiffness matrix is 
obtaine d. 


Compute diaqonal stiffness matrix: 

K ii = M ii j 2Tlt N i ) 2 ' ^fault) 


Read-in full stiffness matrix. 


Integer indicating whether to employ p plane 
approximations to aerodynamics . Functions 
used for this fit have the form: 


2 *j (». + 21. 

5,(p) = A * A P + A 2 p + £ ,, h; 

: 3=1 


J 


1 




Do not compute or use p-plane fit. (Default) 


Employ p-plane curve fit, using coefficients 
which have been previously computed and 
stored on TAP El , 


Compute and employ p-plane approximations for 
full set of aerodynamics from TAPES, 
rewinding TAPE 5 automatically. 


Compute and employ p-plane approximations for 
selected set of aerodynamics from TAPE1 . 


Compute and store geometric coefficients, used 
when computing sensor deflections for plate- 
type sections . 


Do not compute. 


Individual parameters are optional if default values are acceptable. 






tome lis t 


INPUT 


When required 

Aiqehrai c 
variable 

FORTRAN variable 

Descri pt.i on 



I SENSE 




= 1 

Compute and store sensor deflections on TAP El . 



= 0 

Do not compute. 



ICONS YS 

Inte.jer indicating whether control system is to 
be included in analyses. 



■= 1 

Set up control system coefficients (if 
1C ASK - 2 or I SELECT > 0) and perform 
analyses* using a control system. 



--= 0 

Do not i nc lude control system. No quantities 
pertaining to control system need be input. 
( IV* fa ill t ) 



IFLUT 




= 1 

Perform oharact er i st i c root determi nat. i on. 



- 0 

Ho not determine roots for current case. 
(Def aul t) 



I AEROWR 

Inteqer indicatinq whether to write out aerody- 
namic forces to file TAP Hf>. 



- 1 

Write out original aerodynamic force;; if 

I CASK = ±2 and/or selected aerodynamic forces 
if I SELECT / 0. 



- | 

!)o not write out aerodynamic forces. (Default) 



IMODKWR 

Hi i s pa r ame t.e r is e mp 1 oyed primarily in 

non i ni t i a 1 i /.at i on case to output array data 
from previous case. 



- 1 

Write out qenera 1 ized masses, structur.il dump- 
ings, and natural frequencies (generalized 
stiffness if ISTIFF = 1) to file TAP Eh . 



- 0 

Do not write these arrays on TAPK6 unless 
I SELECT * 0 or ICASK = +2. (Default) 



IAP1.T 

Inteqer indicatinq whether to plot aerodynamic 
forces and corresponding curve fits. Uses 

selected model aerodynamics (those on TAPE)). 



- -1 

(Mot only tabular points and interpolated curve 
(refer to KFIT). 



■- 1 

Plot tabular points, interpolated curves 
(refer to KFIT), and p -plane curves (if 
IS PLANE * 0). 



0 

Do not plot aerodynamic data. {Default ) 



IPKPLT 

Inteqer indicatinq whether to generate plots of 
characteristic root data, which type of plots 
to qenerate, and whether to allow modi f i ca- 
t i on if plot parameters during generation of 
plots on a plot -per -plot basis. If 
IPKPLT > 0, modification will be allowed. If 
IPKPLT < 0, only the last set of data will be 
plotted and no changes to plot parameters 
input in namelist. INPUT will be allowed 
during generation of plots. 



tl 

twMier.it e plots of loci versus velocity, alti- 
tude, or density, variation determined by 
KVAR - 1, 2, or 3, respectively, during root 

generation, or versus correspondi ng dynamic 
pressure it KVAR < 0 or versus gain it 
I c ;A I N 1 during root determination. 


t v.i lues .ir e accept able. 


Individual parameters are optional it def aul 


TABLE II 


Set Namelist When required 


Algebraic 

variable 



ISPLANE t 0 



FORTRAN variabl 


Generate plots of C and w versus velocity, 
altitude, or density depending upon KVAR 
during root determination (or corresponding 
dynamic pressure if KVAR < 0) or versus gai 
if IGAIN = 1 during root determination. 

Generate both (C ) and root loci plots, 
n 

Do not qenerate any characteristic root 
plots. (Default) 


Delete unwanted characteristic root data used 
for plotting from data storage file TAP El . 

Do not delete any data. 


Number of coefficients desired currently for 
p-plane fit, including rational terms 
p/ (p + b^J. Maximum = 10. 

Number of polynomial coefficients in p-plane 
fit of jth column of aerodynamic forces. 
(Default = 3) 

Constant denominator coefficient in Ath 

p/(p + . ) corresponding to jth column mode 

If same values of b^^ are used for more than 
one mode, they stillneed only be input once 
See definition of NCOL(J). 

Number of reduced frequencies to be used for 
jth pressure mode fit. Program uses first 
nkk . of total number n^. If 0, the program 
uses all frequencies for the jth mode. 
(Default = 0) 

Integer indicating whether to use a preceding 


Set each b^_. = b^. (Default = 1) 

Input new coefficients for jth mode. 

Integer indicating whether to include Nth 
constraint (as defined below) when fitting 
the Jth column of aerodynamic forces. (See 
appendix A.l 

Nth constraint not employed. (Default) 


Employ Nth constraint 


ICOF ( 1 , J ) = 1 


Force agreement with tabular value for 


IC0F( 2, J) = 1 


Force slope = 


pT Imj_ Qi j (e ) j 


ICOF ( 3 , J ) = 1 For k^ = 0, Force slope = - pi Re [Q^ 2 ( 0} J/bQ^ 


ICOF ( 4 , J) = 1 


0, Force slope 


ICOF( 5 , J ) = l + 1 Force An = 0. 


Qj(k 0 ) = Qj(k 0 
k n 1 is interpol, 


rpolated value of jth 


aerodynamic force at k Q , using a piecewise 
quadratic fit to Q. (k^ (i = 1, ..., n k ). 


Individual parameters are optional if default values are acceptable 








TABLE ir.- Continued 


lame li St 

When required 

Alqebraic 

variable 

FORTRAN variable 

Description 

INPUT 

XSPLANE * 0 


IERPRT 





= 1 

Write percentage errors between p-plune fit 
and actual data. 




r 0 

Do not write out errors. 


rs PLANE * 0 K 
1 K’nK(2,.J) ^ 1 

i 

o 

n 

THETAN 

Normalizin'! factor employed when apply inq 
constraint number 1 in p-plune fit. 
(Refer to appendix A.) (Default = 1) 


ISPLANE t 0 t, 

rcOE(h,J) = 1 

k 

SPKO 

Value of reduced frequency for which p-plane 
approximation must precisely fit; data. 

(See appendix A, constraint number 6.) 


INTEKI' - 1 or 
I SENSE ■- 1 

i 

IS YM 

1 

Modes symmetric with respect to x-axis. 




4 1 

Modes either arit.i symmet r i c or nonsymme t ri c witl 
respect to x-axis. 




NPMX 

Maximum number of nodes in any one structural 
sect i on. 




NSECTNS 

Number of sections into which structure is 
di vi ded. 




N NO! IKS 

Total number of nodes in all structural sec- 
tions at which modal data are available. 




ISE( 1,1) 

Heqinninq node number of structural section I 




1 ss ( I ) 

Undine node number of structural section 1 




I PLATE ( I ) 
- 1 

Indicates that structural section I is modeled 
.is plate. 




. 0 

Indicates that structural section I is modeled 
as beam. 


i 


NPOF 

Number of deqrees of freedom in TAB array. 

If all sections are plate type, this can be 
set to 3 and TAR(I,4) array can be omitted. 
(Default. - 4) 



X 

rq 

XCd 

x-coordi nate of aircraft center of mass. 


! 

X ° 

X0( I ) 

x-eoordinate of translated axis of section I 
when I PLATE ( I ) - 0; x-coordi nat e of mean 

of root chord ot structural section I when 
I PLATE ( I ) 1. 



Y 0 

Y0( I ) 

y-coordinate of translated axis ot section I 
when I PLATE = 0; y-coordinate at root ot 
section 1 for IPLATE(I) - 1. 


i 

. .. i 

/l 

R0( I ) 

Sweep anqle of elastic axis of section I it 
I PLATE ( I ) - 0; if I PLATE ( I ) 1, K0(I) 1/b 

where b R is root, semichord for section I. 


K'UNSYS s 1 

! 

ICS El 'Al> 

Tnteqer indicatinq from where 1 o qet sensor 
deflections needed in conjunction with 
control system. 




1 

Sensor deflections input by user on TAPES. 




l) 

•'“■ilsor deflections .m’ rend-in t r om data 

storaqe tile TAPE 1 (refer to did i nit ion ot 
I SENSE). (Default) 



TABLE II.- Continued 


Name 1 is t 

When required 

Alqebraic 

variable 

FORTRAN variable 

Description 

f INPUT 

ICONSYS = 1 


ICSACT 





-- 1 

Actuator transfer function relates actual con- 
trol surface position to commanded position. 
(Default ) 




= 0 

Cbntrol surfaces are treated as full degrees of 
freedom and actuator transfer functions can 
be considered to relate actuator-produced 
hinge moments to actuator inputs. 



G u 

GN{ I , J) 

Gain for Ith control and Jth sensor pair. The 
user will find it convenient (virtually 
essential in plotting gain root loci) to non- 
dimensional ize . by including desired 

dimensional gain in (T L )^j. Then gain vari- 
ations will be in terms or a multiple of 
a nominal or desired value. (Default = 0) 




PHASE 

Phase error (degrees) in transfer function. 
(Default = 0) 




I PH & JPH 

Control { I PH) and sensor (JPH) pair for which 
phase error is introduced. 




ISCHWTL( I, J) 

1 




= 0 

No scheduling in control logic for control I, 
sensor J. (Default) 




= 1 

Scheduling is included in control logic for 
control I, sensor J. 


ICONSYS = 1 & 
ICASE = + 1 


ICHSEN 
= 1 

Indicates a change is to he made in a denomi- 
nator coefficient for sensor dynamics. (See 

namelist SENINP. ) 




= 0 

Sensor dynamics are unchanged from preceding 
case. (Default) 




ICHFIL 





* 0 

Indicates a change is to be made in filter 
dynamics (^e namelist FILTIN). 




= 0 

No change is to be made in filter dynamics. 
(Default) 




= 1 

Change only occurs in numerator coefficients. 




= 2 

Change only occurs in denominator coef f icients . 




= 3 

Change occurs in both numerator and denominator. 




ICHACT 





* 0 

Indicates a change is to be made in actuator 
dynamics (see namelist ACTINP) . 




= 0 

No change is to be made in actuator dynamics. 




= 1 

Change only occurs in numerator coefficients. 




= 2 

Change only occurs in denominator coefficients. 




= 3 

Change occurs in both numerator and 
denomi nator. 


IFLUT = 1 


KVAR 

Inteqer indicating type of variation to be 
performed in flutter analysis. If negative, 
any corresponding plot output will be with 
respect to dynamic pressure. Can be chanqed 
at plot time in interactive execution. 


individual parameters are optional if default values are acceptable. 



TABLE II. - 


Name list 


When required 


Alqebrai c 
variable 


FORTRAN 


t input 


IFLUT = 1 


KVAR 


variable 


Description 


Velocity is varied with respect to fixed den- 
sity, usinq fixed Mach number aerodynamics. 
RHOO and VO are required inputs. 

Altitude is varied, matching density and 

velocity at the input Mach number. (HO and 
XMACH are necessary inputs, as are the unit 
conversion factors defined subsequently. 


IGAIN 


= 0 


I DMULT 


Density is varied with respect to fixed 

velocity and Mach number. (VO and RHOO are 
necessary inputs.) This option is usually 
run to correlate with model testinq in a wind 
tunnel when only density is varied. 


IV r for in an automatic qain variation. 

(IGONSYS - 1) DV- is set to 0. Velocity and 

density are fixed at VO and RHOO if KVAR = ; 
or +3. Velocity and density are computed t( 
correspond to altitude HO if KVAR = +2. 

Input also DELGAIN. 


Do not perform automatic qain variation. 
( Default ) 


- 1 Multiply all matrix elements by control-ays Lei 

transfer -function common denominator (this 
must be done to avoid sinqu lari ties when 
trackinq control system toots beqinninq wit! 
zero gains) 

-■ 0 Do not clear control denominator. (Default) 


NV 


Number of stepwise variations in quantity beii 
varied - velocity, altitude, density, or 
qain. (Default = 1) 


NIT 


NC1JT 


n 


AV 

c 


N FINE 


Gr iter ion for determininq whether converqence 
upon sufficiently accurate estimate of 
characteristic root has been achieved. If 
s. is Nth estimate of ith root and s. i 
1 N 1 N-1 

(N-1 estimate, then either |s. | < 0.1 f. 

Is. - s, l/ Is. I < e implies conver- 
l L N-1 1/ [ 1 N- 1 \ 

qence. (Default = 10 10 ) 

Maximum number of iterations, without conver- 
qence, to determine roots usinq determinant 
iteration techniques before usinq matrix 
iteration techniques instead. (Default = 3 

Maximum number of times to halve step size wh 
computed step size causes an i ncrease in si 
of determinant. (Default. = 5) 

Number of subdivisions of DV for finer scan. 
This is used when converqence difficulties 
arise in obtaininq a root usinq both deter- 
minant iteration and matrix iteration durir 
normal step size scan. For chanqe in sta- 
bility, a finer scan is also performed. 
Number of subdivisions for this is 
mi n{max(NFINE, 1 0) , 2 5} . (Default - 25) 


I DAMP 


Inteqer which indicates how damping is to 
incorporated into equations of motion 
(eq . ( 1 ) ) . 


Dampinq is J - 1 q s . (Default) 
i 

Damping is (s/|s|)g s (s * 0) and 
ti s (s = o) . 1 


Individual parameters are optional if default values are acceptable. 


be 
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Set Namelist When required 


Alqebraic 

variable 


3 T INPUT IFLUT = 1 


Damping is (s/u) a dq g (if a> n * 0) 

and 0 (if u) n = 0), corresponds to viscous 
damping 2C • 1 

Integer indicating whether to trace Ith 

root, provided Ith initial guess in nonzero. 

Trace root. Use linear (1), quadratic (2), or 
cubic (3) curve fit to extrapolate for esti- 
mate of Ith root at new value of independent 
variable. (Default = 1) 

Do not trace root. 


Print characteristic root data to screen (file 
OUTPUT) every N times they are computed. 
(Default = 0) 


IFLUT = 1 & 

KVAR = +1, +3 


Initial velocity for velocity variation; 
constant velocity for density variation. 


Initial density for density variation; constant 
density for velocity variation. 


IFLUT = 1 & 
KVAR = +2 


Initial altitude for altitude variation. 

Factor for converting input velocity to feet 
per second for use in subroutine AT62. 
(Default = 0.0833 to convert inches per 
second ) 

Factor for converting input altitude to feet 
for use in subroutine AT62. (Default = 1) 


IFLUT = 1 & 
TGAIN = 0 


Factor for converting density, output by sub- 
routine AT62 in slugs per foot 3 , into type of 
units corresponding to input model. 

{Default = 0.48225E-04 for converting to inch 
units) 

Step size for each increment of independent 
variable: velocity, altitude, or density, 

cor responding to KVAR = +1, +2, or +3, 
respective ly . 


IFLUT = 1 & 
IGAIN = 1 


DELGAIN ( I , J ) Increment in GN(I,J) for proper root estima- 

tion. Each nonzero DELGAIN should be equal. 
(Default = 0) 


SII(I) & S2I(I) Two initial estimates (complex) for Ith charac- 
teristic root. S2I(I) should be best guess. 
In some cases, these can be determined auto- 
matically by program. See definition of 
IOPTl and I0PT2. 


Compute initial estimates for elastic modes 
using natural frequencies and FS and GS 
(defined subsequently). (See eq. (17).) 


Do not compute initial estimates in this 


Use matrix iteration to obtain initial estimate 
(appendix D, O ~ 2 and 3). This method 
attempts to converge on elastic and rigid- 
body characteristic roots ( 0 p = 2) unless 
some S2I(I) * 0. (0 p =3), in which case 
method will attempt to converge on nonzero 
roots nearest those input in array S2I(I). 

See also definition of OMR. 


03 









Set Namelist when required 

3 ^ INPUT IPKPLT = ±: 


red 

Algebraic 

variable 

FORTRAN variable 

De: 


FORG, FEND, 
FDEL, & NDF 

Lower limit, upper limit, increment between 
major tic marks, and number of minor divi- 
sions between major tic marks for ordinate 
in plots. (Optional) 

GORD, GEND, 
GDEL, & NDG 

Lower limit, upper limit, increment between 
major tic marks, and number of minor divi- 
sions between major tic marks for ordinate 
in C plots. 

XS(I> & 
YS { I ) 

Arrays defining sensor locations. These should 
be unrotated coordinates I = 1,NS. 

ITYPE ( I) 

Integer indicating whether sensors measure 
linear or angular motion. 

= 1 

Linea r. 

= 2 

Angular. 

NSS( I) 

Number of the structural section on which sensoi 
is located . 




ASD( I, J) 

Ith coefficient of denominator polynomial of 
transfer function which models sensor 
dynamics for Jth sensor (constant term 
first). 1=1, . . . , 11. 

XKS( I, J) 

Ith coefficient of numerator polynomial of 
transfer function which models sensor 
dynamics for Jth sensor (constant term 
first). I = 1 , . . . , 10. 

I ORD ( J ) 

Integer indicating type of sensor that sensor 
is . 

- 0 

Position sensor. (Default) 

= 1 

Rate sensor. 

- 2 

Acceleration sensor. 

TsDYN( J) 

Integer indicating whether to include sensor 
dynami cs . 

= 1 

Include sensor dynamics. 

= 0 

Perfect sensor; do not include dynamics. 


ICASE = ±2 & 
ICONSYS = 1 or 
ICASE = ±1 & 
ICHFIL = 1 


ICASE = +2 & 
ICONSYS = 1 or 
ICASE = ±1 s 
ICHFIL = 1 


K1BD0 
K1BD1 
K1RD2 
AN1 & AN2 


Inteqers indicating control-sensor pair for 
which sensor data are to be filtered. 

I and J refer to identification prior to an 
selection. Only control-sensor pairs for 
which data are to be filtered need be input 
End all filter input with 0,0. 

Natural frequency and damping ratio in second 
order notch. See section "Compensation 
options . " 

Gain for integral feedback. 

Gain for proportional feedback. 

Gain for derivative feedback. 

Time constants' of zero and pole in lead-lag 
filter. 


Individual parameters are optional if default values are acceptabl 










TABLE II.- Continued 


Namelist 

When required 

Alqehrai c 
variable 

FORTRAN variable 

Description 

f F I uTIN 

ICASE = ±2 & 

ICONS YS - 1 or 
ICASE = +1 & 
ICHFIL = 1 

a Fn & a Fd 

A FN ( K , L ) & 
AFD ( K, Ij ) 

Kth coefficient (constant term first) of 
numerator and denominator polynomials for 
Lth filter in transfer function relatinq 
control T to sensor .7. (See set. 6) 

K m 1 , . . . , G. 




r K TLTER 

Array of i nteqers indicat.inq type of filters 
beinq combined to make up transfer function 
relatinq control T to sensor J. Example: 

I FILTER - 3, 4, 1 means transfer function 

will be composed of a type 3 and 4 filter. 
(Defaul t = 1 , no f i 1 ter ) 




= 1 

No further filterinq. (Default) 




= 2 

Notch filter. 




- 3 

Tnteqral filter. 




*. 4 

Proportional plus derivative filter. 




= 5 

T.oad-laq or laq-lead filter. 




--- 6 

General rational function (quotient of two 
polynomi a 1 s ) . 

f ACT INI’ 

ICASE = f2 K 

ICONS YS - 1 or 
ICASE - +1 f. 
1CHACT i 1 

a An * a Ad 

AACTN ( I , J ) K 
AACTP( I, .1) 

Ith coefficient (constant term first) ot 
numeral oi and denominator polynomials of 
transfer function which models actuator for 
Jth cont rol . 




1 ACT(.J) 

Inteqer indicat.inq whether to include actuator 
dynamics* J corresponds to identification 
prior to any selection. 


j 


- 1 

Include actuator dynamics for dth control. 




= 0 

Perfect actuator; do not include any dynamics. 
( Defaul t ) 

* SELECT 

I SELECT * 0 


NMODES 

Total number of modes to bn selected from 
available data. After selection, NM is set 
to NMODES. 




NRNEW 

New number of riqid-body modes; i.e., new NR. 




NCNEW 

New number of control modes; i.e., now NC. 




NIX' 

Asoendinq array of mode numbers indieatinq for 
which modes modal data are to he selected. 
(NMODES = Total number of nonzero NOC values 
Example: NOC = 2, 7, and 9 and NMODES - 3 

indicates that data for three modes, name ly 
nodes 2, 7, and 9, will be selected for 

analysis. This option effects a truncation 
of the mode l . 




NKNEW 

New number of reduced frequencies; i.e., new 
NK. 




= 0 

Indicates no frequency selection; all will he 
used. (Default) 


I SELECT * 0 6, 
NKNEW / 0 


NOK 

Ascend inq array of reduced frequency indices 
indicat.inq for which frequencies nerodynnmii 
force data are to be selected. Example: 

NOK =1, 2, 5 and NKNEW = 3 indicate only 

aerodynami c force data cor respondi nq to thu 
reduced frequencies, namely 1st, 2nd, and 
5th, will be used in analyses. 


individual parameters are optional it default values are acceptable. 




r 


TABLE II.- Continued 


Set 

Namelist 

When required 

Algebraic 

variable 

FORTRAN variable 

Description 



IDELPLT * 0 


IPLT2 

Array of numbers indicating which sets of 

characteristic root data are to be saved from 
file. If ICASE > 0 (interactive execution) 
then these can be modified before actual 
deletions occur. 

1 1 


IPKPLT > 0 


ICOUNT 

Number of sets of characteristic root data 
to be plotted, currently, out of total 
number of sets of data generated during 
different cases. Complete list available is 
printed out before this input is required. 

■ 

HI 

IPKPLT > 0 


IPLT 

Array of numbers indicatinq which sets of 
characteristic root data are to be plotted. 

1 

■ 

IPKPLT > 0 


ICOMBIN 

Inteqer indicating whether to combine plots of 
different sets of characteris tic root data 
into one plot (all on same frame) or generate 
plots on separate frames. 

I 




= 1 

Combine plots. 

■ 

1 



= 0 

Do not combine plots. 

1 4 


IPKPLT > 0 & 
ICOMBIN = 1 


PTITLE 

80-character description of data being plotted. 
This will appear on plot being generated in 
lieu of oriqinal plot titles when plots are 
being combined. 

1 5 

+ PLOTP 

IPKPLT > 0 


I DASH 

Indicates type of line used to connect points 
in data currently being plotted. Ibis may be 
changed for each set of data plotted. 





= 0 

(solid line) 





= 1 
- 2 






= 3 






= 4 






= 5 






= 6 






= 7 






ISAME 

Integer indicating how many of next sets of 
characteristic root data can be plotted with 
same PLOTP namelist data currently being 
read-in, including current set. The program 
will not request another PLOTP namelist to be 
input until n - ISAME sets of data have been 
plotted. (Default = 1) 





KVAR 

Input negative of value of KVAR that was used 
when determining roots if you desire to 
change the independent variable in the 
graphics output (see definition of KVAR on 

p. 101). 





MODEPK( I) 

Integer indicating whether to generate plot 
showing Ith characteristic root variation. 





= 1 

Generate plot of data. (Default) 





= 0 

Bypass Ith root information. 



IPKPLT = +1 or +3 


XRLORG, XRLEND, 
XRLDEL, & 
NDXRL 

Lower limit, upper limit, increment between 
major tic marks, and number of minor divi- 
sions between major tic marks for abscissa 
in root-locus plots. (Optional) 


individual parameters are optional if default values are acceptable. 
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TABLE II.- Concluded 


Name list 

When required 

Alqebrai c 
variable 

FORTRAN variable 

Description 

+ PLOTP 

IPKPLT m +1 or +3 


YRLORG, YRLEND, 
YRLDFJL, & 
NDYRL 

Lower limit, upper limit, increment between 
major tic marks, and number of minor divi- 
sions between major tic marks for ordinate 
in root -locus plots. (Optional) 


IPKPLT = +2 or +3 


VORD, VEND, 
VDEL, & NDV 

Lower limit, upper limit, increment between 
major tic marks, and number of minor divi- 
sions between major tic marks for abscissa 

in C and w plots. (OptionaL) 

n 




FORG, FEND, 
FDEL, &, NDF 

Lower limit, upper limit., increment between 
major tic marks, and number of minor divi- 
sions between major tic marks for ordinate 
in plots. (Optional) 




GORG, GEND, 

GO EL, & NOG 

Lower limit, upper limit, increment between 
major tic marks, and number of minor divi- 
sions between major tic marks for ordinate 
in C plots. 


IPKPLT > 0 


IOP 

Inteqer indi ca tinq whether to continue 
plotting. This option allows user to 
qene rate additional plots of characteristi 
root data that either were not included 
before or for which changes in the plots a 
des i. red. 




= 1 

Continue plot tinq. 




-■= 0 

Return to main program and continue with a n 
case or terminate execution. 


f Individual parameters are optional if default values are acceptable. 




TABLE III.- DESCRIPTION OF TAPE5 INPUT DATA 


[All data are read-in using free-field format] 




Algebraic 

variable 

FORTRAN 
array name 

Description 

READ statements 

1 



KK 

Reduced frequencies at which 
aerodynamic data were 
computed 

FOR K=1 ,NK 
READ* , KK( K) 
where KK(K)=k v 

Xv 

2 

ICASE = ±2 
and/or 
I SELECT > 0 

Re {Q } & 
Im{Q} 

VARDR & VARDI 
or 

AR & AI 

Real and imaginary parts of 
aerodynamic motion forces; 
pairs are read-in bv columns 
(per pressure mode ) f per 
frequency and transposed 
internally 

FOR K=1 , NK 

FOR J=1 , NM ( +NG ) + 

FOR 1=1 ,NM 

READ*, ( VARDR (K,J, I) , VARDI (K,J, I ) 
where VARDR ( K , J , I ) =Re { Q 0c„ ) _ T } 

tv 1 r J 

and VARDI (K, J,I)=Im{p(k ) } 

K I , J 

3 

ICASE = ±2 
and/or 
ISELECT > 0 

M 

MAT1 

Generalized masses; read-in 
by rows 

FOR 1=1 ,NM 
FOR J=1 , NM 
READ* , MAT1 (I, J) 
where MAT1 (I, J)=M I 

1 

ICASE = ±2 
and/or 
ISELECT > 0 

< f n> 

FREQ 

Natural frequencies (Hz) 

FOR 1=1 , NM 
READ*, FREQ(I) 
where FREQ(I)=f 

n I 

5 

ICASE = ±2 
and/or 
ISELECT > 0 
& 

ISTIFF - 1 

K 

MOMSQ 

Generalized stiffness matrix; 
read-in by rows 

FOR 1=1 ,NM 
FOR J=1 , NM 
READ*, MOMSQ ( I, J) 
where MOMSQ( I , J) =KI , J 

6 

ICASE = ±2 
and/or 
ISELECT > 0 

V 

G 

Structural damping coefficients 

FOR 1=1, NM 
READ* , G ( I ) 
where G( I ) =g 

S I 

1 

ICASE = ±2 
and/or 
ISELECT > 0 
& 

ICS READ = 1 

H 

DS 

Sensor deflections; read-in by 
columns per sensor 

FOR J=1 ,NS 

FOR 1=1 , NM 

READ*, DS ( I, J) 

where DS(I,J)=H^ T 
I , J 

8 

ICASE = ±2 
and/or 
ISELECT > 0 
& 

ICSACT = 0 


HE 

Actuator elongation per unit 
generalized coordinate; 
read-in by rows per control 

FOR 1=1, NC 
FOR J=1 , NM 
READ* , HE( I , J) 

9 

ICASE = ±2 
and/or 
ISELECT > 0 
& 

ICSACT = 0 

f d 

FD 

Generalized force per unit 

actuator hinge moment; read-in 
by columns per control 

FOR J= 1 , NC 
FOR 1=1 , NM 
READ * , FD ( I , J ) 


+ If gust forces are included as last NG columns of aerodynamic data, provision has been made to 
read these in (NG^O), although they are not stored in these arrays. 
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TABLE IV." DESCRIPTION OF TAP El 0 INPUT DATA 


[All data are read-in using free-field format] 


Set 

When 

required 

A1 gebraic 
variable 

FORTRAN 
array name 

Description 

READ statements 

1 

INTERP = 1 or 
ISENSE = 1 

x or y' 

TAB( I , 1 ) 

x-coordinates if I PLATE = 1 
or points along y'-axis 
if I PLATE = 0 at which 
modal data are available 

FOR 1= 1 , NNODES 
READ* , TAB ( I , 1 ) 
where TAB(I,1)=Xj or y^. 

2 

INTERP = 1 or 
ISENSE = 1 

y or z (y ' ) 

TAB (1,2) 

y-coordinates at which 

modal data are available 
if I PLATE = 1 or bending 
deflections along the 
y'-axis if IPLATE = 0 

FOR 1=1, NNODES 
READ* , TAB (1,2) 
where TAB{I,2)=y I or z(y^.) 

3 

INTERP = 1 or 
ISENSE = 1 

z or Ry (y ’ ) 

TAB ( I , 3) 

Modal deflections at (x, y) 
if IPLATE = 1 or rotation 
about y' at points along 
y'-axis if IPLATE = 0 

FOR 1=1, NNODES 
READ*, TAB( I, 3) 
where TAB(I,3)=z I or R^(yJ-) 

1 

INTERP = 1 or 
ISENSE = 1 
and 

NDOF = 4 

0 or R^(y' ) 

TAB ( 1 , 4) 

Zeros if IPLATE = 1, or 
rotation about x' at 
points along y'-axis if 
IPLATE = 0 

FOR 1=1, NNODES 
READ* , TAB( 1 , 4) 
where TAB(I,4)=0 or R^(y' ) 


Repeat sets 1 to 4 for each mode (1 to 3 if NDOF = 3) 
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TABLE V.- DATA STORED ON DATA STORAGE FILE TAPE1 BY PROGRAM 


Record 

Algebraic variable 

FORTRAN array name 

Size 

Description 

1 

M 

MAT1 

NMODES* *2 

Selected generalized masses 

2 


MOMSQ 

NMODES 

2 

Selected Mu)^ array 

3 

f n 

FREQ 

NMODES 

Selected frequency (Hz) array 

4 

S [bjy] 

SCOF & BN 

NCOEF*NMODES* *2 + 250 

Selected p-plane coefficients and 
b^j array 

5 


IA1SUB 

102 

Subindex array for PKPLOT data 

6 

{k} T & Re [Q(k i ) f i = 1 , . . . ,n k ] 

RFRQ & AR 

NK+NK*NMODES**2 

Selected reduced frequencies and 
real part of aerodynamic data 

7 

Im[Q(k i ) , i=1 , • . • ,n k ] 

AI 

NK*NMODES* *2 

Selected imaginary part of 
aerodynamic data 

8 

»E 

HE 

NCNEW* NMODES 

Selected actuator elongation per 
unit generalized coordinate 

9 

f d 

FD 

NMODES *NCNEW 

Selected generalized force per 
unit actuator hinge moment 

10 

H 

CS 

NS* (NMODES -NCNEW ) 

Selected sensor deflections 

1 1 


G 

NMODES 

Selected structural dampings 

1 2 


DCT, CNTLN, & ICLN 

NCTM+ ( MPM+ 1 ) *NS*NCNEW 

Control logic common denominator 
(d(s)), numerator elements, 
and number of coefficients in 
each numerator element 

14 

{%j> s ib ij 5 

mam 

NC0EF*NM* *2+250 

Full set of p-plane coefficients 
and b^j array 

15 

! 

INDEX 4 for SPLINE 

51 

Subindex array for geometric 
coefficients used in surface 
spline interpolation 

16 

H 

CS 

NS* (NM-NC ) 

Full set of sensor deflections 

1 7 


INDEX for FILTIN 

51 

Subindex used in storing filter 
data 

20 

a Sn' a Fn' 

ASN, ACLN, . . . 

Computed 

All intermediate control data 
used in computing d and N^ 

50 


All labeled commons 


All labeled common data from 
previous case for restart 
capability 
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Subroutine SYKSU°r 


Compute geometric 
ddtii for symmetric 


cose. 















2 


Jf I S PLANE - 3 » call SPLANE 
overlay to determine p- 
plane fit to aerodynami cs 
for selected modes. 


If I CONS Y S = 1 and 
!CASE= - 2 or i f I CONSYS- 1 
and some changes are to 
be made to control -system 
dynamics nr if 1 CONS YS = I 
and ISELLCTfSO, CONTROL 
overlay is called. In 
CONTROL overlay matrix 
of transfer functions 
relating sensor input to 
actuator output is 
determined . 




Program SPLANE 


Sets up array 
storage areas. 



Enter \ 
AERCOEF / 




Subrout i r 

e AERCOEF 

Selected aerodynai 
in fit are read-i 
constrained least- 
performed to detet 

Coefficients are ; 

lie data employed 
from TAPE 1 . A 
-squares fit. is 
-mine coefficients. 

itored on “APE 1 . 

l 


{ Ij.'-OL ) 

Program CONTROL 

Determine sensor 
dynamics . 


Determine filter 
dynamics . 


Determine actuator 
dynamics . 


Overal 1 matrix of 
transfer functions 
is obtained by 
performing matrix 
product of sensor 
f i 1 ter and 
actuator transfer 
function matrices . 

These data are 
stored on random 
access file TAP E 1 
for use in sub- 
sequent overlay. 



.Subroutine SENOYN 


Compute sensor 
transfer function 
Ji 

times s w for 
each sensor (see 
name! ist SEN1NP) . 


Enter 

CLOGIC 


D 


[Subroutine CLOGIC | 

Polynomial 
representat ion 
| o f control logic | 

transfer- 

functions . 


^ / Enter \ 
! FILDYN / 


Subr 

outi ne FI 

LDYN 

Cent! 

rol 1 

o n i c 


tran: 

s fer 

func 

tions 

are 1 

bui 1 1 

up 

by 

input set 

ect i 

O n 

from 

bui 1 

t - i n 


filter types . 


(See 

namel i st 


tilt 

IN.) 




< Enter' \ 
ACTBYN / 




Subroutine ACTDYN 

Actuator transfer 
functions are 
obtained using 
the data from 
name list ACTINP. 

' 



Figure 5.- 


Continued • 
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Outboard aileron 


Vertical accelerometers 



Figure 6.- Paneling of wing and horizontal tail for aerodynamic force computation 


Filter 


(i ; 2 - e 

|, g units 

i' 2 2 2 2 

s J 11 250 s + 40s + 196“ s + 220s + 456 

6 l , dec, 

1 

s + 2s + D s s + 250 ^2 + 1QOs + 190 _; 3 2 s 2 + 400s + 52Q 2 



-1 i G 1X = (3.147 a 10 J ) (l.61 ' (q - 1034) < 1 


75 < D = 735 - 833. 3N + 0. 0033267a < 300 
~ S Ma 


if N < 0.61, set N = 0.61 
Ma Ma 


') y 

if q v 23.924 kN/itT, set q = 23.9248 kN/m 


deg 


Actuator 


14 

6.583 10 


/ 2 



6 \ 

/ 2 

_ „ , 6 \ 

I s + 681.7 

s + : 

l .02* 

* * L0 ) 

f s + 654. 

6s + 1.597 ' 10 j 


vS, 


deq 


Figure 7.- Control law and actuator dynamics. 




Imaginary 







Re ( \ j f rad/sec 


(a) altitude root loci. 

Figure 9.- Open-loop stability characteristic variation with 

altitude at = 0.86. 








Re( ) , rad/sec 

(a) Density root loci. 

Figure 10.- Open-loop stability characteristic variation with density. 

N = 0.86; U = 269.6 m/sec. 











Velocity, m/sec 


(a) Open loop. 

Figure 12.- Characteristic root variation with velocity (reduced-order model) 

N =0.86; P = 0.6101 kg/m 3 . 











Figure 14.- Effect of phase errors in 
loci at U ( p-p analysis) 



, j , rad/sec 


ntrol 


signal elastic mode gain root 
0.86; h = 4.572 km. 







Locus of short-period root with altitude 
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